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Introduction 

fl'/j 

The  >rk  described  here  represents  a  combined  theoretical  and  experimental  effort 

to  address  the  mechanical  behavior  of  interfaces  in  brittle  fiber-reinforced  brittle  matrices. 

In  the  three  year  program,  it  was  the  aim  of  the  study  to  examine  variables  which  allow  for 
toughness  optimization  of  ceramic  fiber-reinforced  ceramic  matrix  composites.  The  results 
highlighted  below  range  the  gambit  from  the  development  of  analytical  solutions  to  solve 
the  eigenstrain  problem  of  both  a  single  fiber  (treated  as  a  cylindrical  inclusion)  in  an  elastic 
half-space  to  the  experimental  verification  of  toughening  models  for  one  fiber-reinforced 
system.  We  have  further  examined  the  validity  of  a  variety  a  test  techniques  which  are  used 
to  acquire  the  material  properties  for  composite  systems,  and  have  designed  a  new  single 
fiber  pullout  technique  reminiscent  of  actual  fiacture  in  fiber-reinforced  composites. 

Finally,  additional  tests  on  model  fiber  systems  have  been  performed  to  measure 
parameters,  such  as  bridging  forces,  necessary  for  toughening  models. 

II.  Highlights 

(A)  Inclusion  Method  Applied  to  Fiber-Reinforced  Composite  Problems 

Exact  solutions  have  been  developed  for  axisymmetric  stresses  and  displacement 
fields  caused  by  a  solid  cylindrical  inclusion  in  the  presence  of  a  uniform  eigenstrain  in  a 
half  space  by  using  Green's  functions.  Green's  functions  were  obtained  for  problems  of 
an  elastic  half  space  with  a  free  surface  or  rigidly  fixed  surface  for  a  homogeneous  isotropic 
elastic  solid  in  earlier  work.[  1-3]  The  strain  energy  is  also  obtained  in  closed  form.  These 
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solutions  are  now  of  the  form  which  allows  for  evaluation  of  load  diffusion  and  absorption 
capabilities  in  the  matrix  as  a  function  of  fiber  length  and  material  properties.(See  Appendix 
6.)  Such  solutions  are  necessary  to  examine  problems  of  fiber  debonding  and  pullout  from 
a  matrix. 

(B)  Theoretical  Fracture  Maps 

An  energy  method  has  been  used  to  examine  the  amount  of  toughening  which 
results  from  the  bridging,  debonding,  and  sliding  processes  in  ceramic  matrix  ceramic  fiber 
composites.  Using  the  equivalent  inclusion  method,  the  crack  is  approximated  as  a  thin 
spheroid  bridged  by  a  series  of  fibers.  In  prior  studies,  this  method  was  used  to  calculate  a 
"fracture  map"  which  described  conditions  for  crack  propagation  and  arrest  with  and 
without  fiber  reinforcement.  The  maps  demonstrate  the  various  contributions  of  bridging, 
debonding  and  frictional  sliding  for  both  constant  and  Coulomb  friction.  Most  recently, 
this  method  has  been  extended  to  examine  debond  initiation  stresses  and  bridging  stresses 
as  a  function  of  interface  toughness,  residual  thermal  stresses,  volume  fraction  and  fiber 
radius.(See  Appendix  3.)  These  stresses  are  compared  to  the  fiber  strength  and  hence, 
describe  conditions  where  debonding  and  pullout  are  favored  over  fiber  fracture. 
Calculations  of  bridging  stresses  are  compared  with  those  measured  in  the  experimental 
program  described  below. 

(C)  Bridging  Forces  by  Inclined  Fibers  or  Whiskers 

In  general,  fracture  toughening  models  neglect  any  fibers  which  are  not  aligned 
parallel  to  the  axis  of  applied  stress  under  the  supposition  that  they  contribute  little  to  the 
toughening  process.  In  order  to  test  this  assumption,  model  composite  specimens 
containing  fibers  of  different  inclination  angle  were  fabricated  with  SiC  monofilaments  and 
borosilicate  glass  matrices  via  hot-pressing.  The  double  cleavage  drilled  compression 
fracture  mechanics  specimen  was  used  to  characterize  crack  bridging  by  the  inclined  fibers. 
During  stable  crack  growth,  the  crack  opening  displacement  is  monitored  using  a  computer 
data  acquisition  system.  The  bridging  force  was  extracted  using  a  fracture  mechanics 
analysis  involving  finite  element  computation.  Results  are  shown  in  Appendix  1  and  these 
represent  the  first  quantitative  experimental  measurements  of  crack  bridging  by  inclined 
fibers  in  ceramic  composites.  At  small  crack  opening  displacements,  the  bridging  forces 
decrease  with  inclination  angle.  However,  the  bridging  forces  remain  significant  for  large 
inclination  angles  over  long  distances  (many  orders  of  magnitude  of  the  fiber  diameter) 


indicating  that  the  fibers,  though  undergoing  a  bending  stress,  do  not  automatically  fracture 
in  the  crack  wake  region  .  This  is  most  likely  due  to  matrix  spalling,  shown  in  Appendix 
1,  which  diminishes  the  severity  of  the  bending  angle.  These  results  are  useful  for 
toughening  models  of  both  fiber-  and  whisker-reinforced  systems  where  the  fiber 
reinforcement  is  not  perpendicular  to  the  crack  plane. 

(D)  Fiber  Debonding  and  Pullout  using  a  Novel  Pullout  Test 

In  order  to  develop  a  basic  understanding  of  the  fiber  debonding  and  pullout 
process  as  a  function  of  interfacial  chemistry  and  interfacial  roughness,  a  number  of 
interfacial  (pullout)  tests  were  evaluated.  None  were  satisfactory  to  to  adequately  represent 
crack  growth  through  the  matrix  of  a  real  composite  coupled  with  fiber  interception, 
debonding,  and  pullout.  In  order  to  correct  this  deficiency,  a  modified  single  fiber  pullout 
test  was  designed.  In  this  test,  a  stiff  fixture  holds  the  specimen  in  place:  a  single  fiber 
embedded  along  its  entire  length  in  matrix  material.  A  pre-crack  is  started  in  the  matrix  and 
grown  to  intercept  the  fiber.  The  resultant  load  displacement  is  shown  in  Appendix  2.  The 
curve  shows  the  corresponding  non-linear  load  displacement  curve  is  indicative  of 
progressive  debonding  of  the  fiber,  followed  by  complete  debonding  and  finally,  frictional 
sliding.  Pullout  lengths  of  over  10  mm  have  been  observed  in  the  SiC 
monofilament/borosilicate  glass  system.  Using  the  analytical  methods  of  Kerans  and 
Parthasarathy  [4]  and  Marshall  [5],  an  interfacial  fracture  toughness  (Gi)  of  approximately 
0.2  J/m^  and  a  coefficient  of  friction  (p)  of  0.14  have  been  measured. 

(E)  Uniaxial  Tensile  Tests  of  Fiber  Composites 

In  order  to  compare  experimental  results  with  earlier  crack  arrest  theory  (Section 
B),  studies  on  fiber-reinforced  cements  have  been  conducted.  Specifically,  the  stress- 
strain  response  of  steel  fiber-reinforced  cements  with  varying  fractions  of  fiber- 
reinforcement  were  measured.  Of  particular  interest  is  the  composite  stress  at  the  bend- 
over-point  in  the  stress-strain  curve.  The  theoretical  treatment  establishes  a  critical  volume 
fraction  of  reinforcement  below  which  debonding  precedes  the  bend-over-point  in  the 
stress-strain  response  and  above  which  nonlinearity  occurs  prior  to  debonding.  Good 
agreement  is  demonstrated  between  experiment  and  theory. 
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The  main  points  from  this  work  which  could  be  used  by  the  scientific  and 
engineering  communities  are  as  follows: 

1 .  Solutions  were  developed  for  stress  and  displacement  fields  for  circular  cylindrical 
inclusions  adjoining  the  surface  of  a  half-space.  These  solutions  could  be  extended 
to  treat  single  fiber  pullout  from  a  matrix. 

2.  A  new  single  fiber  test  technique  has  been  designed  to  determine  interfacial 
properties  of  fiber-reinforced  materials  that  more  accurately  reflects  behavior  in 
full-scale  composites  than  any  previous  single-fiber  test. 

3.  The  first  quantitative  measurements  of  bridging  forces  were  collected  for  inclined 
fibers  in  ceramic  composites.  Such  measurements  are  imperative  for  accurately 
modelling  toughening  in  fiber-reinforced  and  whisker-reinforced  materials. 

4.  "Fracture  maps"  were  designed  to  establish  conditions  for  crack  arrest,  debonding, 
pullout,  etc.  in  continuous  fiber-reinforced  ceramics. 
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Abstract 

Crack  bridging  by  inclined  fibers  has  been  studied  in  a  brittle  fiber-brittle  matrix 
model  ceramic  composite.  Results  of  the  fiber  bridging  force  vs.  the  crack  opening 
displacement  have  been  obtained  for  different  fiber  inclination  angles  using  a  fracture 
mechanics  approach.  Localized  matrix  cracking  at  the  fiber  has  been  observed  and 
related  to  fiber  inclination.  The  experimental  results  showing  the  effect  of  the  fiber 
inclination  angle  are  discussed  and  compared  with  theoretical  analyses  to  provide 
insight  into  crack  bridging  by  inclined  fibers/whiskers.  Implications  for  toughening  by 
whisker  bridging  will  also  be  discussed.  [Key  words:  crack  bridging,  fibers,  whiskers, 
inclination,  composites] 
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L  Introduction 


Significant  increases  of  fracture  toughness  can  be  achieved  by  incorporating 
whiskers  into  ceramic  matrices.^'^  Studies  of  whisker-reinforced  ceramic  composites 
have  established  that  the  major  toughening  mechanism  in  these  composites  is  crack 
bridging  by  whiskers.®'®  A  study  of  large-grain  polycrystalline  aluminas  shows  a 
similar  toughening  mechanism:  crack  bridging  by  grains.®*^®  These  studies^'^®  have 
analyzed  the  mechanics  of  crack  bridging  by  whiskers  and  grains.  Due  to  the  lack  of 
information  regarding  the  dependence  of  the  bridging  force  on  the  whiskert^  inclination 
angle,  the  reinforcements  are  assumed  to  be  perpendicular  to  the  crack  plane^’®’^®  or 
a  simplistic  relation  between  the  bridging  force  and  the  inclination  angle  is  used.^^ 
Alternatively,  the  inclined  whiskers  are  neglected  because  they  are  considered  to  fail 
by  bending  and  result  in  negligible  toughening.®  In  reality,  the  whiskers  can  be  of  any 
orientation  with  respect  to  the  crack  plane  in  whisker-reinforced  composites.  The 
large  msyority  of  whiskers  are  not  perpendicular  to  the  crack  plane.  A  more  realistic 
model  would  consider  whiskers  of  all  orientations.  Indeed,  recent  in  situ  observations 
for  a  SiC-whisker  reinforced  alumina  have  indicated  the  importance  of  inclined 
whiskers  to  the  bridging  process.  Therefore,  it  is  essential  to  explore  the  effects  of 
the  fiber/whisker  inclination  angle  on  the  degree  of  crack  bridging  by  each  individual 
fiber/whisker. 

There  have  been  a  number  of  studies  on  the  effect  of  fiber  inclination  angle  on  fiber 
pullout  in  systems  other  than  brittle  fiber-brittle  matrix  composites.^®'^^  However, 
pullout  of  inclined  fibers  in  brittle  fiber-brittle  matrix  composites  has  received  very 
limited  treatment  until  recently.  Wetherhold^^  analyzed  the  pullout  of  a  short  brittle 
fiber  poorly  bonded  to  the  matrix.  The  simple  analysis  indicates  that  the  frictional 
ptillout  force  is  inversely  proportional  to  the  cosine  of  the  inclination  angle  0  of  the 
fiber  with  respect  to  the  siuface  normal  for  a  given  fiber  embedded  length.  This 
relationship  is  assumed  to  apply  for  0  less  than  a  physically  plausible  critical  angle 
beyond  which  the  fiber  will  fi:acture  due  to  bending  stresses. 

Recently  Leung  and  Li^®  conducted  an  extensive  finite  element  analysis  to  study 
the  effect  of  fiber  inclination  angle  on  the  crack  bridging  stress  in  brittle  fiber-brittle 
matrix  composites.  Matrix  spalling  was  incorporated  into  the  analysis  based  upon 


^'hrhe  term  whiskers  may  also  refer  to  elongated  grains  and  short  fibers. 


2 


observations  in  ducllie  fiber-brittle  matrix  composites.  The  results  for  long  fibers 
show  that  the  bridging  stress  decreases  with  an  increase  of  the  inclination  angle  for  a 
given  crack  opening  displacement  (COD).  The  maximum  crack  opening  displacement 
befo:  fiber  failure  is  smaller  for  inclined  fibers  than  aligned  fibers.  Nevertheless,  the 
bridg^g  stress  firom  inclined  fibers  remains  significant,  particiilarly  in  the  near  crack- 
tip  region  where  the  crack  opening  displacements  are  small. 

In  spite  of  this  progress  in  the  theoretical  aspect,  there  is  little  experimental  work 
on  the  effect  of  fiber  inclination  on  the  fiber  pullout  or  crack  bridging.  The  objective  of 
this  work  is  to  characterize  experimentally  crack  bridging  by  inclined  fibers/whiskers. 
Based  upon  the  results  for  a  single  inclined  whisker,  the  toughening  contribution  due 
to  crack  bridging  of  fibers/whiskers  of  distributed  orientations  can  be  estimated.  The 
results  are  expected  to  be  useful  in  providing  insight  into  the  mechanism  of  whisker 
bridging  and  in  improving  the  present  models  of  whisker  bridging. 


IL  Fracture  Mechanics  Approach 


For  the  present  work,  we  chose  to  examine  the  effect  of  fiber  inclination  angle  by 
studying  a  fiber  bridging  a  crack  which  closely  simulates  the  conditions  of  crack 
bridging  in  an  actual  composite.  The  double  cleavage  drilled  compression  (DCDC) 
fracture  mechanics  specimen  was  used  for  this  purpose. 

The  DCDC  specimen  was  developed  by  Janssen  for  studying  slow  crack  growth  in 
glasses.  The  applied  stress  intensity  Kapp  for  the  DCDC  specimen  is  related  to  the 
measured  crack  length  (a)  and  the  applied  compressive  stress  ((T^)  as  follows: 


O’, 

*PP  "  g(o/r) 


(1) 


where  r  is  the  central  hole  radius,  g(a/r)  is  the  dimensionless  geometric  fiuiction  of  the 
crack  length  (a)  and  specimen  geometry.  The  geometric  function  can  be  determined 
experimentally  by  testing  specimens  of  known  toughness  or  obtained  by  theoretical 
analyses.  The  advantages  of  the  DCDC  configuration  (Fig.  1)  include  simple 
specimen  shape  and  loading  geometry,  stable  crack  propagation  over  a  large  region  of 
specimen  length,  and  the  fact  that  there  is  no  need  fw  crack  guiding  gro 
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The  DCDC  specimen  has  been  used  to  study  the  crack  interaction  with  an  aligned 
fiber^®  and  to  quantify  the  toughening  effect  of  fibers  with  a  ductile  interfacial  layer 
The  possibilities  for  quantifying  the  effects  of  bridging  fibers  using  the  DCDC 
specimen  were  demonstrated  in  terms  of  the  traction  exerted  across  the  crack. As 
shown  schematically  in  Fig.  2(A),  the  crack  is  bridged  by  a  fiber  behind  the  crack  tip. 
The  influence  of  bridging  can  be  expressed  in  terms  of  the  change  in  the  stress 
intensity: 


A/fb  ~  ■  -^app  (2) 

where  Kc  is  the  base  toughness  of  the  matrix  and  /fapp  is  the  applied  stress  intensity 
(Negative  AA'b  indicates  a  closing  action  on  the  crack;  positive  A^b  indicates  an 
opening  action  on  the  crack).  Such  results  are  useful  for  examining  the  toiighening 
effect  of  a  single  fiber  or  whisker.  However,  for  quantitative  analyses  of  the 
propagation  of  a  bridged  crack,  the  key  information  is  the  relation  between  the 
bridging  force  and  the  crack  opening  displacement. 

The  crack  opening  displacement  (COD)  can  be  measured  directly  using  a 
displacement  gauge.  The  bridging  force  must  be  extracted  from  the  measured  change 
in  the  stress  intensity.  It  is  necessary  to  establish  the  relation  between  the  change  in 
the  stress  intensity  and  the  bridging  force  in  the  DCDC  fracture  mechanics  specimen. 
The  bridging  action  of  a  fiber  behind  the  crack  tip  can  be  approximated  by  self- 
equilibrated  traction  forces  {Pit)  acting  at  a  bridging  distance  (6)  behind  the  crack  tip, 
as  illustrated  in  Fig.  2(B).  For  bridging  forces  very  close  to  the  crack  tip  (small 
bridging  distance)  the  near  tip  solution  can  be  used.  For  general  values  of  the  bridging 
distance,  the  crack  bridging  relation  in  the  DCDC  specimen  is  not  available.  A  finite 
element  analysis  was  used  to  obtain  the  relationship  between  the  change  in  the 
stress  intensity  and  traction  forces.  In  this  analysis,  the  change  in  the  stress 
intensity  (A/0  was  computed  as  a  fimction  of  the  bridging  force  {Plt\  the  crack  length 
(a),  and  the  bridging  distance  (6).  For  the  given  DCDC  specimen  geometry,  the 
normalized  change  in  stress  intensity  [A/fVrAP/O]  is  a  function  of  the  normalized 
crack  length  {dr)  and  the  normalized  bridging  distance  (6/r).  Typical  resiilts  obtained 
by  the  finite  element  analysis  are  shown  in  Fig.  3  for  three  values  of  normalized  crack 
length.  From  these  numerical  calculations,  the  crack  bridging  force  was  computed 
firom  the  experimentally  measured  change  in  the  stress  intensity. 
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m.  Experimental  Procedure 


The  model  composite  system  consists  of  SCS-6^  silicon  carbide  fibers  and  a 
borosilicate  glass^^  matrix.  Specimen  fabrication  procedures  were  similar  to  those 
used  by  Butler  et  al.^  The  fibers  were  degreased  in  1,1,1-tncbloroethane  and  then 
rinsed  with  acetone  and  ethanol.  The  glass  plates  were  washed  with  acetone  and 
ethanol.  Both  the  fibers  and  the  glass  plates  were  then  dried  at  about  llO^C  for  60 
minutes.  The  fibers  were  sandwiched  between  two  glass  plates.  The  sandwich 
assembly  was  placed  in  a  hot  press, $  which  was  then  evacuated  using  a  mechanical 
pump.  The  temperature  was  raised  to  730^^0  before  the  hot  pressing  pressure  was 
increased  slowly  to  1  MPa  and  held  for  60  minutes.  Specimens  with  fiber  inclination 
angles  from  0  to  40  degrees  were  fabricated.  Each  specimen  with  two  fibers  was 
machined  to  a  DCDG  specimen  with  the  fibers  located  symmetrically  on  both  sides  of 
the  center  hole.  There  were  four  specimens  for  each  fiber  inclination  angle.  The 
specifications  of  the  DCDC  specimen  dimensions  are  given  in  Fig.  1.  For  present 
work,  the  DCDC  specimen  with  the  Janssen  geometry  {HIW  =  10,  W/r  =  3.75)^®  was 
used  with  the  center  hole  radius  of  0.795  mm  and  the  specimen  thickness  of  3.5  mm. 

Before  testing,  each  DCDC  specimen  was  precracked  at  the  top  and  the  bottom 
edges  of  the  hole  where  the  cracks  would  nucleate  and  grow.S^  A  Knoop  indenter  with 
an  indentation  load  of  29.4  N  was  used  for  this  puri)ose.  Compression  testing  of  the 
DCDC  specimen  was  carried  out  on  an  Instron  universal  testing  machine^  at  a 
crosshead  speed  of  0.025  mm/min.  The  testing  was  conducted  in  a  controlled 
environment  of  flowing  dry  nitrogen  gas  to  minimize  the  effect  of  slow  crack  growth  on 
the  measured  crack  resistance  of  the  glass.  A  crack  opening  displacement  (COD) 
gauge  was  attached  to  measure  the  COD  at  the  fiber  location.  A  computerized  data 
acquisition  system  was  used  to  collect  the  crack  length,  load  and  crack  opening 
displacement  data.  The  crack  length  was  monitored  using  a  travelling  microscope 
connected  to  a  digital  micrometer.  Once  the  crack  tip  was  located,  the  crack  length, 
load  and  COD  data  were  collected  simultaneously.  A  data  point  was  collected  every 

^Silicon  Carbide/Carbon  fiber,  AVCO  Specialty  Materials,  Lowell,  MA. 

**Soda  borosilicate  glass.  Coming  Code  7740,  Coming  Glass  Works,  Coming,  NY. 

^Vacuum  hot  (.''ess,  GCAA^acuum  Industries,  Inc. 

S^It  is  not  necessary  to  pre-crack  the  DCDC  tincture  mechanics  ^specimen.  To  the  authors' 
knowledge,  none  of  previous  studies  pre-cracked  the  DCDC  specimen.  However,  pre-cracking  has  the 
advantage  of  ensuring  that  cracks  to  starting  at  the  center  line.  Without  pre-cracking,  the  crack  can 
propagate  unstably  to  much  greater  distances  than  with  pre-cracking,  particularly  if  tJie  load  train 
is  not  very  stiff. 

Ilnstron  Model  1125,  Instron,  Corp.,  Canton,  MA. 
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few  seconds.  Tests  were  stopped  when  the  crack  grew  to  a  length  of  13  to  15  times  of 
the  center  hole  radius.  The  maximum  COD  at  the  fiber  was  in  the  order  of  10  pm. 
Borosilicate  glass  specimens  without  fibers  were  tested  to  obtain  the  dimensionless 
geometric  fimction  g(a/r)  for  the  DCDC  specimen  and  for  comparison  purpose. 

After  mechanical  testing,  the  specimen  was  bonded  with  a  cyanoacrylate-based 
adhesive^^  to  a  glass  slide  and  a  section  of  the  specimen  containing  one  of  the  fibers 
was  cut  with  a  slow  speed  diamond  saw.  After  the  adhesive  was  removed  with 
acetone  in  an  ultrasonic  bath,  two  pieces  of  matrix  glass  held  together  by  a  bridging 
fiber  were  obtained.  These  two  pieces  of  glass  were  then  separated  for  microscopic 
examination.  The  fracture  svirface  near  the  fiber  was  examined  using  scanning 
electron  microscopy  (SEM)  to  observe  possible  matrix  spalling  related  to  pullout  of 
inclined  fibers. 


IV.  Results 

The  experimentally  determined  dimensionless  geometric  function  of  the  DCDC 
specimen  is  shown  in  Fig.  4  for  the  crack  length  range  of  interest.  Results  were 
obtained  from  two  tests  and  provide  almost  identical  geometric  functions.  This 
calibration  fiuiction  was  used  to  compute  the  applied  stress  intensity  for  all  other 
tests  using  Eq.  (1). 

The  applied  stress  intensity  vs.  the  crack  length  is  presented  in  Fig.  5  for  the 
cases:  (A)  without  fibers,  (B)  fibers  at  (p  =  0®,  (C)  fibers  at  0  =  24®,  (D)  fibers  at 
0  =  40®.  In  the  case  without  fibers,  the  applied  stress  intensity  is  nearly  independent 
of  the  crack  length  and  is  equal  to  the  base  toiighness  of  the  borosilicate  glass.  In  the 
presence  of  bridging  fibers,  the  applied  stress  intensity  increases  as  the  crack 
intercepts  and  passes  the  fiber,  as  shown  in  Fig.  5(B),  (C)  and  (D).  The  results 
indicate  a  trend  of  decreasing  bridging  force  with  increasing  inclination  angle.  Another 
interesting  observation  is  that  the  applied  stress  intensity  tends  to  decrease  slightly 
as  the  crack  approaches  and  passes  the  fiber  before  ifapp  increases  again  for  the 
case  of  large  inclination  angles.  This  was  consistently  observed  in  a  number  of  tests 
with  large  fiber  inclination  angles  and  is  not  well  imderstood.  Fiber  inclination  might 
have  induced  a  mode  II  stress  intensity  which  promotes  crack  growth. 

^iDuro™  Quick  Gel™  super  glue,  Loctite  Corp.,  Cleveland,  OH. 
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Based  upon  the  applied  stress  intensity  data,  the  equivalent  bridging  force  (as 
defined  in  Fig.  2(B))  was  computed  using  the  finite  element  analysis  described  briefly 
in  Section  II.  The  crack  bridging  information  consists  of  the  bridging  force  as  a 
fimction  of  the  corresponding  crack  opening  displacement.  Examples  of  the  bridging 
force  vs.  COD  are  presented  in  Fig.  6  for  fiber  inclination  angles  of  0,  8,  24,  and  40°. 
For  small  fiber  inclination  angles,  the  bridging  force  levels  off  at  a  COD  value  of 
approximately  8  pm,  indicating  that  fiber  firictional  sliding  over  the  embedded  length 
starts  at  this  point.  Comparison  of  these  results  shows  a  trend  of  decreasing  bridging 
force  with  an  increase  in  the  fiber  inclination  angle  in  the  range  of  COD  studied,  as 
expected  from  the  stress  intensity  results. 

The  bridging  force  is  plotted  as  a  function  of  the  inclination  angle  for  different 
crack  opening  displacements  in  Fig.  7.  The  results  show  a  general  trend  that  the 
bridging  force  decreases  with  increasing  fiber  inclination  angle  for  crack  opening 
displacements  between  2  and  10  pm.  No  data  is  available  for  greater  crack  opening 
displacements  due  to  two  limitations:  1)  specimen  size  and  2)  the  crack  length  at 
which  the  test  was  stopped  to  prevent  specimen  shattering  at  high  compressive 
stresses. 

SEM  micrographs  of  the  matrix  at  the  fiber  are  shown  in  Fig.  8  for  fiber  inclination 
angles  of  8°,  24°  and  40°.  For  the  fiber  inclination  angle  of  8°,  fiber  pullout  occurred  as 
the  two  pieces  of  matrix  glass  were  separated.  Little  matrix  cracking  was  observed 
(Fig.  8(A)).  For  large  inclination  angles,  the  fiber  failed  as  the  two  pieces  of  matrix 
glass  (held  together  by  the  fiber)  were  separated  in  the  direction  normal  to  the  crack 
plane.  Localized  matrix  cracking  was  observed  (Fig.  8  (B)  and  (O).  The  observed 
cracking  pattern  looks  somewhat  similar  to  cracks  generated  by  Vickers  indentation. 
The  severity  of  matrix  cracking  increases  with  the  fiber  inclination  angle,  although 
part  of  the  matrix  cracking  apparently  has  resulted  from  SEM  sample  preparation. 
In  another  observation,  the  matrix  glass  pieces  were  separated  in  the  fiber  direction 
using  a  fixture  in  order  to  minimize  matrix  cracking  due  to  SEM  sample  preparation. 
This  observation  revealed  much  less  severe  matrix  cracking.  Further,  the  lateral 
crack  did  not  intercepted  the  surface. 
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V.  Discuission 


We  have  studied  the  effect  of  fiber  inclination  on  crack  bridging  in  a  SiC  fiber¬ 
glass  matrix  model  composite  system.  The  results  show  that  the  bridging  force  by 
inclined  fibers  decreases  with  increasing  fiber  inclination  angle  in  the  range  of  crack 
opening  displacement  observed.  This  finding  is  in  qualitative  agreement  with  Leung 
and  Li's  theoretical  analysis.^^.  Leung  and  Li  have  shown  that  the  total  crack 
bridging  force  consists  of  contributions  due  to:  1)  the  bending  force  on  the  inclined  fiber 
and  2)  fiber  pullout.  The  crack  bridging  force  is  the  sum  of  these  two  contributions 
in  a  direction  normal  to  the  crack  plane.  The  contribution  from  fiber  bending 
increases  with  the  fiber  inclination  angle  initially  but  decreases  at  higher  angle.  The 
effect  of  localized  matrix  cracking  relax  the  more  heavily  stressed  fibers  and  thus 
reduce  the  difference  between  fibers  lying  at  different  angles.  The  fiber  pullout 
contribution  in  the  direction  normal  to  the  crack  plane  decreases  with  the  inclination 
angle.  Leung  and  Li^^  pointed  out  that  for  most  composites  the  bridging  force  is 
usuaUy  dominated  by  the  fiber  pullout  component.  Therefore,  the  total  bridging  force 
decreases  with  the  fiber  inclination  angle  (see  Ref  18  for  details). 

It  is  expected  that  bending  of  inclined  fibers  eventually  leads  to  fiber  failure  as  the 
crack  opening  displacement  increases.  Nevertheless,  the  present  study  shows  that 
the  bridging  force  for  inclined  fibers  is  nowhere  near  zero.  For  example,  crack  bridging 
by  fibers  inclined  by  40°  is  sustained  more  than  6  mm  behind  the  crack  tip  (or  more 
than  40  times  the  fiber  diameter).  Further  crack  propagation  and  large  COD's  would 
be  required  for  fiber  fracture.  In  actual  whisker-reinforced  materials,  the  majority  of 
the  whiskers  are  inclined  with  respect  to  the  crack  plane.  Therefore,  bridging  by 
inclined  whiskers  is  expected  to  make  up  the  main  portion  of  toughening  in  the  neeu* 
tip  wake  region.  If  the  present  results  can  be  extended  to  whisker-reinforced  ceramic 
composites,  toughening  predictions  based  on  aligned  fibers  would  represent 
overestimates.  In  addition,  the  dependence  of  bridging  force  on  the  whisker  inclination 
angle  and  crack  opening  displacement  also  has  a  profound  effect  on  the  the  shape  of 
the  crack  resistance  of  the  composite.  Compared  with  predictions  with  aligned 
whiskers,  the  slope  of  the  initial  part  of  crack  resistance  curve  would  be  smaller  if  the 
bridging  force  at  a  given  COD  is  reduced  as  a  result  of  whisker  inclination. 

Comparison  of  matrix  cracking  between  smaller  and  large  fiber  inclination  angles 
in  Fig.  8  shows  that  localized  matrix  cracking  is  related  to  fiber  inclination.  It  is 
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observed  that  larger  fiber  inclination  angles  result  in  a  greater  degree  of  matrix 
cracking.  Matrix  cracking  reduces  the  severity  of  fiber  bending  and  has  the  benefit  of 
delaying  fiber  breakage.^*  Accordingly,  the  fiber  can  maintain  integrity  a  greater 
distance  behind  the  crack  tip  as  a  result  of  matrix  cracking.  At  the  same,  matrix 
cracking  also  modifies  the  force-displacement  relation. 

The  present  observation  that  the  crack  bridging  force  decreases  with  the  fiber 
inclination  appears  to  be  in  contradiction  with  Wetherhold's  simple  analysis^  ^  which 
concludes  that  the  bridging  force  should  increase  with  the  fiber  inclination  angle.  It 
must  be  pointed  out  that  the  results  of  the  present  study  are  obtained  for  the  case 
where  the  fiber  pullout  dominates  the  total  bridging  force.  Such  cases  include  those 
with  large  embedded  fiber  length  (e.g.,  the  embedded  length  more  than  20  times  of  the 
fiber  radius  for  present  study)  and  high  interfacial  shear  strength  and  firictional  shear 
stress.  For  very  short  whiskers  or  elongated  grains,  or  very  low  interfacial  shear 
strength  and  frictional  shear  stress,  increased  firictional  shear  or  interlocking  due  to 
whisker  inclination  could  dominate  the  bridging  force.  In  this  case,  the  bridging  force 
is  expected  to  increase  as  a  results  of  whisker  inclination  angle  until  whiskers  fail  or 
the  matrix  spalls.  Indeed,  high  values  of  the  coefficient  of  firiction  (greater  than  unity) 
required  for  parametric  fitting  of  crack  resistance  curves  in  some  alumina  ceramics 
using  a  grain-bridging  model  can  be  attributed  to  geometrical  irregularities  including 
grain  interlocking.^^  Based  upon  the  foregoing  discussion,  it  is  clear  that  there  is  no 
general  trend  for  the  effect  of  the  fiber/whisker  inclination  angle  on  crack  bridging.  It 
is  determined  by  other  factors  in  addition  to  the  fiber/whisker  inclination  angle.  The 
present  study  covers  a  limited  range  of  these  factors,  and  further  studies  are  required 
to  obtain  a  complete  quantitative  picture. 
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Fig.  1.  Schematic  of  the  DCDC  specimen  in  which  2H  is  the  specimen  height,  2W 
is  the  specimen  width,  t  is  the  thickness,  r  is  the  radius  of  the  center  hole,  a 
is  the  crack  length,  and  is  the  applied  compression  stress.  For  the 

Janssen  geometry,  HlVf  =  10  and  W/r  =  3.75. 

Fig.  2.  Schematic  of  crack  bridging  in  the  DCDC  specimen:  (A)  Fibers  at  an 
inclination  angle  0  bridge  the  DCDC  cracks;  (B)  Representation  of  the 
bridging  force  by  traction  force  Pit.  The  distance  from  the  crack  tip  to  the 
bridging  force,  or  the  bridging  distance,  is  denoted  h. 

Fig.  3.  Crack  bridging  relations  in  the  DCDC  fracture  mechanics  specimen 
showing  the  normalized  change  in  the  stress  intensity  as  a  fimction  of  the 
normalized  bridging  distance  for  different  crack  lengths.  The  results  were 
obtained  using  a  finite  element  analysis.  The  dimensions  of  the  DCDC 
specimen  are:  HI  W  =  10,  W/r  =  3.75,  and  r  =  0.795. 

Fig.  4.  The  dimensionless  geometric  function  of  the  DCDC  fracture  specimen 
obtained  experimentally. 

Fig.  5.  Examples  of  the  results  showing  the  applied  stress  intensity  vs.  the 
normalized  crack  length  for  the  cases:  (A)  without  fibers,  (B)  fibers  at 
0  =  0®,  (C)  fibers  at  0  =  24®  and  (D)  fibers  at  0  =  40®. 

Fig.  6.  Examples  of  the  results  showing  the  bridging  force  vs.  the  crack  opening 
displacement  for  different  fiber  inclination  an^es. 

Fig.  7.  The  bridging  force  as  a  fimction  of  the  inclination  angle  for  different  crack 
opening  displacements.  Open  symbols  indicate  data  points  from  a  single 
test. 

Fig.  8.  Matrix  cracking  related  to  fiber  inclination  at  inclination  angles:  (A)  8®,  (B) 
24®,  and  (040®. 
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FIBER  DEBONDING  AND  PULLOUT  PROCESSES 
IN  CERAMIC  COMPOSITES 
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Evanston,  IL  60208-3108 

The  fiber  debonding  and  pullout  behavior  in  a  model  ceramic  composite  system  has 
been  investigated  using  a  novel  single-fiber  pullout  technique.  Stable,  progressive 
debonding  was  observed  prior  to  fiber  fracture  and  frictional  fiber  pullout.  Interfa¬ 
cial  parameters  were  calculated  from  single  load- displacement  curves  using  a  recent 
analysis  of  the  fiber  debonding  and  pullout  process.  The  interfacial  parameters 
extracted  from  these  experiments  are  consistent  with  independent  measurements 
reported  in  the  literature.  The  morphology  of  the  fiber  surface  appears  to  play  a 
significant  role  in  determining  the  fiber  debonding  and  sliding  behavior. 

INTRODUCTION 

It  is  now  acknowledged  that  the  mechanical  properties  of  fiber-reinforced  ceramic 
composites  are  controlled  primarily  by  the  interfacial  properties.^’^  Recognition  of 
the  importance  of  the  interfacial  properties  has  prompted  the  development  of  a 
number  of  techniques  to  evaluate  these  parameters,  including  the  fiber  pullout^’* 
and  pushdown’  tests.  The  work  described  here  utilizes  a  variant  of  the  fiber  puUout 
technique. 

Extraction  of  interfacial  properties  from  single-fiber  pullout  tests  requires  fitting 
of  the  experimental  data  to  appropriate  micro-mechanics  based  models  of  the  fiber 
debonding  and  pullout  process.  A  number  of  such  analytical  models  have  been 
developed  recently  with  varying  assumptions  and  approximations.  A  recent  analysis 
by  Kerans  and  Parthasarathy,^  which  has  been  used  to  analyze  fiber  push-down 
data,*  is  also  applicable  for  the  fiber  pullout  experiment.  The  model  is  based  on 
progressive  debonding  of  the  fiber  and  accounts  for  the  thermal  mismatch  of  the 
fiber  and  matrix,  including  the  residual  axial  stress,  and  specifically  addresses  the 
role  of  fiber  surface  asperities. 

The  model^  considers  a  fiber  of  radius  r  embedded  in  a  cylinder  of  infinite  radius 
of  the  matrix  with  a  length  lo  protruding  from  the  matrix.  The  response  load  Pa  is 


related  to  tlie  displacement  at  the  end  of  the  fiber  during  progressive  debonding  as 
follows: 


IpPa  1  -  2i/fk 
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where  k  is  given  as: 
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and  Ef  and  Em  ue  the  elastic  moduli  of  the  fiber  and  matrix,  Vf  and  Vm  are  the 
respective  Poisson’s  ratios  and  \i  is  the  friction  coefficient.  The  parameter  Pr  is  the 
residual  axial  load  on  the  fiber  due  to  thermal  expansion  mismatch  between  the 
fiber  and  matrix: 


=  (3) 

where  Aa  is  the  thermal  expansion  coefficient  difference  for  the  matrix  and  fiber 
(om  -  OE/)  and  AT  is  the  change  in  temperature  over  which  the  residual  stresses 
develop. 

The  term  Pa  corresponds  to  the  critical  axial  load  in  the  fiber  at  the  debond  crack 
tip  necessary  to  propagate  the  debond  crack.  It  is  directly  related  to  the  critical 
str^  energy  release  rate  of  the  interface  according  to: 


'  ~  4ir»r3£, 

P*  is  a  critical  axial  load  on  the  fiber  at  which  Poisson’s  contraction  completely 
counteracts  the  residual  normal  stress  and  the  surface  asperity  effects  such  that 
there  is  no  frictional  resistance  to  sliding.  When  the  response  load  Pa  equals  P*, 
the  debond  crack  will  propagate  through  the  specimen  at  constant  load  with  a  fixed 
length  frictional  sliding  zone  behind  the  debond  crack  tip.  P*  is  directly  related  to 
the  effective  normal  stress  at  the  interface,  ,  as  follows: 


P*  = 


(5) 


It  is  our  intent  in  this  study  to  use  a  unique  single-fiber  pullout  technique  to 
experimentally  investigate  the  progressive  debonding  behavior  in  a  model  ceramic 
composite.  The  role  of  the  interfacial  friction,  interfacial  debond  energy,  residual 
stress  state  and  fiber  surface  morphology  will  be  explored.  The  results  are  to  be 
analyzed  with  the  model  discussed  above  to  extract  interfacial  properties  from  the 
progressive  debonding  portion  of  individual  load-displacements  curves. 


EXPERIMENTAL  PROCEDURES 


For  this  investigation,  single-fiber  composite  specimens  were  fabricated  from  SiC 
monofilaments*  embedded  in  a  soda-borosilicate  glass^  The  SiC  monofilaments  are 
produced  by  chemical  vapor  deposition  (CVD)  of  SiC  onto  a  SS^tm  diameter  carbon 
core.  The  CVD  process  results  in  a  high  degree  of  surface  roughness,  which  has  been 
shown  to  have  a  profound  effect  on  the  fiber  debonding  and  sliding  behavior.^  The 
surface  of  the  monofilament  consists  of  two  carbon  rich  layers  which  are  believed  to 
provide  the  weak  interface  necessary  for  fiber  debonding  and  pullout. 

The  samples  were  prepared  by  hot-pressing  in  an  argon  atmosphere  at  tempera¬ 
tures  of  approximately  750°C  for  15  minutes  under  a  stress  of  5.0  MPa.  To  fabricate 
samples  with  the  geometry  discussed  below,  it  is  necessary  to  have  the  embedded 
fiber  parallel  to  the  faces  of  the  hot-pressed  disk.  To  accomplish  this,  one  half  of  the 
glass  frit  was  cold  pressed  into  the  hot-press  die  under  a  stress  of  approximately  3.0 
MPa.  The  top  ram  was  removed  and  a  segment  of  the  monofilament  was  placed  into 
the  die.  The  remaining  glass  frit  was  then  cold  pressed  into  the  die  at  an  equivalent 
stress.  The  die  was  then  loaded  into  the  hot-press  and  processed  according  to  the 
conditions  outlined  above. 

The  sample  geometry  that  was  developed  and  adopted  for  use  in  these  experi¬ 
ments  is  shown  in  Figure  1.  A  rectangular  block  of  the  matrix  contains  a  single  fiber 
running  down  the  axis  of  the  sample.  The  sample  was  notched  to  have  a  greatly 
reduced  cross-section  at  the  midpoint  of  the  sample. 

This  modified  single  fiber  pullout  geometry  was  chosen  to  provide  a  more  realistic 
modeling  of  the  fiber  debonding  and  pullout  process  during  composite  fracture  rela¬ 
tive  to  the  standard  single  fiber  pullout  geometry  used  in  previous  studies.^’”  During 
the  testing  procedure,  a  crack  is  propagated  through  the  matrix  leaving  the  fiber  as 
a  bridging  element.  Therefore,  an  artificial  surface  is  not  used  to  model  the  matrix 
crack  and  any  artificial  barrier  to  debond  initiation  is  avoided.  The  interaction  of 
the  matrix  crack  with  the  fiber  to  initiate  debonding  is  more  clearly  represented. 

The  machined  samples  were  directly  bonded  to  the  load  train  of  a  standard 
mechanical  test  framed  using  a  cyanoacrylate  based  adhesive*,  taking  care  not  to 
bond  the  exposed  fiber  ends,  as  indicated  in  Figure  2. 

During  the  pullout  test,  the  matrix  crack  opening  displacement  (COD)  was  mea¬ 
sured  using  a  traveling  microscope  attached  to  a  digital  micrometer.^  The  microm¬ 
eter  was  interfaced  to  a  data  acquisition  system  which  simultaneously  recorded  the 
load.  The  crosshead  speed  was  held  constant  at  lO^m/min. 

*SCS-6  SiC  Monofilament,  Textron  S]>ecialt]r  Materials,  Lowell,  MA 

'Pyrex  #7740,  Coming  Glass  Works,  Coming,  NY 

^Instron  4505  testing  machine,  Instron  Corporation,  Canton,  MA 

iDnro  Quick  Gel,  Loctite  Corporation,  Cleyeland  OH 

*Mitntoyo  di^matic  micrometer  head,  MTI  Corporation,  Paramns,  NJ 


RESULTS  AND  DISCUSSION 


Representative  data  for  these  experiments  are  shown  in  Figure  3,  where  two  sets 
of  data  are  superimposed.  The  load  increases  continuously  with  the  crack  opening 
displacement  until  fiber  fracture.  The  non-linearity  in  the  initial  portion  of  the  data 
is  expected  as  the  compliance  of  the  sample  increases  with  progressive  debonding. 
After  the  sharp  drop  in  load,  the  load  decreases  further  with  decreasing  embedded 
length  as  the  debonded  fiber  is  frictionally  pulled  out  of  the  matrix. 

It  is  interesting  to  note  that,  although  the  peak  loads  and  the  COD’s  at  the  peak 
load  for  the  two  data  sets  are  very  different,  the  final  pullout  lengths  are  nearly  equal. 
Furthermore,  the  magnitude  of  the  peak  loads  and  the  large  load  drop  with  fiber 
fracture  indicate  that  the  debond  crack  had  progressed  far  ahead  of  the  eventual 
fiber  fracture  point. 

Figure  4  demonstrates  the  application  of  the  Kerans  and  Parthasarathy  model^ 
to  the  data,  showing  the  saturation  load  P*  as  well  as  the  debond  propagation  load 
Pd.  Fiber  fracture  occurs  well  before  the  saturation  load  P*  is  reached. 

The  Kerans  and  Parthasarathy  model  (Equation  1)  w‘as  fit  to  individual  data  sets 
using  a  three-parameter  non-linear  least  squaresll  fitting  routine.  For  the  sample 
geometry  used  herein,  the  initial  length  lo  is  zero  and  the  first  term  of  Equation  1 
is  omitted.  The  fitted  parameters  ft,  P4,  and  P*  for  each  of  the  data  sets  shown  in 
Figure  3  are  reported  in  Table  1.  From  these  fitted  parameters,  the  effective  normal 
stress  and  the  interfacial  strain  energy  release  rate  were  calculated,  and  are  also 
reported  in  Table  1. 


Table  1:  Fitted  parameters 


Data  Set 

P<iW 

P*  (N) 

(MPa) 

1 

0.144 

1.17 

104.9 

201.0 

0.24 

2 

0.138 

0.89 

192.4 

0.14 

The  strain  energy  release  rate  is  found  to  be  very  low,  consistent  with  the  weak 
interface  expected  in  these  materials.  The  friction  coefficient  agrees  well  with  the 
measurements  of  Parthasarathy,  et  al.*  for  fiber  pushdown  experiments  in  similar 
materials  and  the  measurement  of  ft  for  graphite  sliding  on  glass  by  Gupta.^ 

The  calculated  effective  normal  stress  accounts  for  both  the  residual  thermal 
mismatch  stress  and  the  fiber  surface  asperity  effects,  according  to  the  following 
relation:^ 


-EmEf 


Ef{l  +  I/,n)  +  Eni(l  —  Vf) 

where  A  is  the  fiber  surface  asperity  amplitude. 


(AaAr+:^)  (6) 

The  radial  thermal  expansion 
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coefficient  for  the  SiC  monofilament  hat  been  measured  by  Goettler  and  Faber^  to 
be  2.63  X  10~*  /"C.  Thking  AT  to  be  the  strain  point  of  the  glass  minus  ambient,  and 
with  a  thermal  expansion  coefficient  for  the  borosilicate  glass  of  3.5  x  10~^  /"C,  the 
residual  normal  stress  for  this  system  is  calculated  from  the  first  term  of  Equation  6 
to  be  20.7  MPa.  The  fiber  surface  asperity  amplitude  can  then  be  determined 
from  Equation  6.  The  asperity  amplitudes  calculated  for  data  sets  1  and  2  are 
0.26  and  0.25|un,  respectively.  These  values  agree  very  well  with  the  independent 
measurements  of  Jero,  et  al.^°  who  used  an  optical  interference  technique  to  directly 
measure  the  surface  roughness.  Their  results  indicate  that  the  average  fiber  surface 
asperity  amplitude  is  0.09-0. 18/xm  with  larger  asperities  up  to  0.36/un. 

The  fiictional  pullout  portion  of  the  load-displacement  curve  can  be  used  to  check 
the  consistency  of  the  fitted  parameters  obtmned  from  the  progressive  debonding 
portion  of  the  data.  The  post-fiber-fracture  embedded  length  is  related  to  /i,  P*, 
and  the  applied  load  by  the  following  relation:^ 


If  the  friction  coefficient  is  taken  to  be  equal  to  that  measured  from  the  progres¬ 
sive  debonding  data,  the  effective  normal  stress  at  the  debonded  sliding  interface 
can  be  calculated.  The  embedded  length  during  pullout  is  approximately  equal 
to  the  final  pullout  length,  measured  after  the  test  is  complete,  minus  the  matrix 
crack  opening  displacement.  The  measured  pullout  lengths  and  the  P*  values  cal¬ 
culated  from  the  data  immediately  following  fiber  fracture  are  reported  in  Table  2. 
The  P*  values  obtained  from  the  parameter  fit  of  the  progressive  debonding  are 
shown  for  comparison.  As  P*  is  a  direct  measure  of  the  effective  normal  stress,  very 
good  agreement  is  found  in  the  interfacial  parameters  obtained  from  the  progressive 
debonding  data  and  the  frictional  pullout  data. 


Table  2:  Post-fiber-fracture  analysis 


Data  Set 

Pullout  Length  (ftm) 

P;ro,  (N) 

Pfric^  (N) 

1 

762 

109.4 

2 

728 

116.5 

HiiBi 

SUMMARY 

A  novel  modified  single-fiber  pullout  technique  was  utilized  to  investigate  the  fiber 
debonding  and  sliding  behavior  in  a  model  ceramic  composite  system.  Progressive 
debonding  was  observed,  consistent  with  recent  analytical  models.  The  morphology 
of  the  fiber  surface  appears  to  have  a  significant  effect  on  the  fiber  sliding  behavior, 
both  during  progressive  debonding  prior  to  fiber  failure  and  during  frictional  fiber 
pullout  fdlowing  fiber  fracture. 


The  Kerans  and  Parthasarathy  analysis^  was  used  to  extract  interfacial  parame¬ 
ters  from  the  progressive  debonding  portion  of  single  load-displacement  curves.  The 
data  from  these  experiments  agree  with  the  model  predictions,  and  the  interfacial 
parameters  calculated  from  the  experimental  data  are  consistent  with  independent 
measurements  reported  in  the  literature. 
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Glass  Matrix 


Figure  1:  Modified  single-fiber  pullout  specimen  geometry. 


Figure  2:  Modified  single-fiber  pullout  test  setup. 
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Figure  3:  Load-displacement  data  for  the  modified  single  fiber  pullout  test. 
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Figure  4:  Application  of  the  Kerans  and  Parthasarathy  model^  to  the  experimental 
load-displacement  data. 
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TENSILE  FRACTURE  OF  FIBER-REINFORCED  COMPOSITES: 
INFLUENCE  OF  THERMAL  RESIDUAL  EFFECT 

W.  B.  Tsai  and  T.  Mura 

Department  of  Civil  Engineering  and  Theoretical  and  Applied  mechanics, 
Northwestern  University,  Evanston,  IL  60208,  U.S.A. 

Abstract-The  inclusion  model  is  used  to  analyze  the  tensile  fracture  of 
fiber-reinforced  brittle-matrix  composites  containing  many 
unidirectionally  aligned  continuous  fibers  with  uniform  strength.  During 
the  debonding  of  fiber-matrix  interfaces,  the  fibers  slide  against  the 
friction  which  is  caused  by  Poisson  effect  and  thermal  expansion 
mismatch.  In  general,  three  fracture  mechanisms  including  interfacial 
debonding,  matrix  cracking  and  fiber  failure  may  occur  in  some  sequence. 
The  sequences  are  strongly  influenced  by  the  thermal  residual  effect  due 
to  temperature  change  AT.  The  analysis  provides  explicit  crack 
length/applied  stress  relation  for  each  mechanism.  Fracture  maps  are 
constructed  to  analyze  the  influence  of  thermal  residual  effect  on  the 
fracture  sequences  and  final  composite  strength.  Two  categories  of 
fracture  transitions  are  defined  analytically  for  the  large  crack  limit. 

1.  INTRODUCTION 

In  a  brittle  matrix  composite  containing  unidirectionally  aligned 
fibers,  three  fracture  mechanisms  are  commonly  expected,  including 
debonding  of  fiber-matrix  interfaces,  matrix  cracking  and  fiber  failure. 
Using  the  standard  fracture  mechanics  approach,  Marshall  et  al.  (1985, 
1987),  Budiansky  et  al.  (1986)  and  McCartney  (1987)  analyzed  the 
influence  of  material  parameters  such  as  frictional  shear  stress,  fiber 
strength  and  interfacial  fracture  toughness  on  the  fracture  mechanisms. 
The  mechanics  of  the  crack  bridged  by  fibers  can  be  also  analyzed  by  the 
inclusion  method  (Mura,  1987).  Mori  and  Mura  (1984,  1988)  investigated 
the  effect  of  fibers  on  the  arrest  of  a  crack  for  the  complete  bridging  and 
complete  sliding  cases.  Yang  et  al.  (1991)  extended  this  analysis  to  the 
case  of  partially  debonded  fiber-matrix  interfaces  assuming  a  constant 
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friction  v/as  present  in  the  debond  region.  Tsai  and  Mura  (1991)  further 
considered  the  origin  of  friction  which  is  the  Poisson  effect  and  thermal 
expansion  mismatch.  Similar  subjects  were  studied  by  Gao  et  al.  (1988), 
Sigl  and  Evans  (1989),  Hutchinson  and  Jensen  (1990).  In  this  paper,  the 
authors  emphasize  the  influence  of  thermal  residual  effect  on  each 
fracture  mechanism.  For  instance,  the  threshold  stress  for  debonding 
initiation  increases  as  the  thermal  residual  effect  is  reduced.  The  reason 
will  be  discussed  later. 

The  present  analysis  also  provides  the  fracture  criterion  and  the 
unique  relation  between  crack  length  and  critical  applied  stress  for  each 
kind  of  fracture  mechanism.  Consequently,  it  is  possible  to  construct  the 
fracture  maps  shown  in  Fig. 8  and  Fig. 9.  Sequences  of  fracture 
mechanisms  and  final  composite  strength  are  easily  obtained  by  drawing 
crack  extension  paths  in  the  fracture  maps.  Also,  the  large  crack  limit  will 
be  discussed  at  the  end  of  this  paper.  Two  categories  of  fracture 
transitions  will  be  defined  and  critical  conditions  are  given  in  terms  of 
material  properties. 


2.  MODEL  AND  ANALYSIS 

2.1  Inclusion  model  and  stress  analysis 

Consider  a  composite  containing  many  aligned  and  randomly 
distributed  continuous  fibers  of  radius  r  and  average  spacing  X. 
Perpendicular  to  the  fiber  direction  is  a  matrix  crack  (Fig.l).  The  shape  of 
the  matrix  crack  is  approximated  by  a  thin  oblate  spheroid  as 

xl  c 

Q-  ■  ^  ,  2.  +  ^<i  ^  ii«l.  (1) 

a 

a  is  the  radius  of  a  penny  shape  crack.  Cq  is  the  half  of  the  fictitious 
thickness  of  the  crack.  The  absolute  value  of  Cg  does  not  appear  in  the 
relevant  equations  to  be  derived  later.  The  cross  section  of  a  fiber  in  the 
crack  is  also  approximated  as 

x'^+xl  xl  c 

i2:  '  ,  ■?  +^<1  ,  -2-«l.  (2) 
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When  an  uniform  tensile  stress  at  infinity  is  applied  in  the  direction 
normal  to  the  crack  surface,  the  eigenstrain  ejj  for  the  equivalent 
inclusion  method  is  introduced  to  simulate  the  bridged  crack.  Specially 


=  ««/> 
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(4) 


Here,  a  is  the  parameter  characterizing  the  frictional  sliding  of  fibers 
along  the  fiber-matrix  interfaces.  There  are  two  extreme  cases 
illustrating  the  physical  meaning  of  a.  First,  complete  bridging  (a  =  0,  no 
sliding  of  fibers)  occurs.  Second,  complete  sliding  of  fibers  (a  =  l)  takes 
place.  The  general  expression  of  a  will  be  derived  later  when  we  discuss 
the  crack  wake  frictional  sliding  of  fibers.  According  to  the  averaging 
processes  employed  by  Mori  et  al.  (1983),  one  obtains  the  average  internal 
stresses  in  fig  -Q,  and  Qs,  given  as 
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Here,  <  >  denotes  the  average  defined  in  the  respective  domain,  p  is  the 
shear  modulus  and  v  is  the  Poisson  ratio  of  the  composite.  /  is  the  volume 
fraction  of  fibers  and  equals  to  .  To  satisfy  the  traction  free  condition 
on  the  crack  surface 

cr  +((T.  \  ,  =0,  (7) 

A  \  33/^2  -Qs 


e„  is  obtained  as 


2(1 


v)a  a 
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The  total  stress  of  the  fiber  in  Q,  (bridged  domain),  sum  of  the  external 
and  internal  stresses,  is  written  as 


K)a; 

Using  Eq.(5),  (8)  and  (9),  the  total  stress  (bridging  stress)  is  given  below: 
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2.2  Analysis  of  frictional  sliding  of  Hbers 
(determination  of  a) 

In  this  section,  the  authors  investigate  the  mechanism  of  frictional 
sliding  of  fibers  to  determine  a.  For  simplicity,  one  considers  a  rod- 
cylinder  model  consisting  of  a  rod  (fiber)  embedded  in  a  cylinder  (matrix), 
Fig.2.  In  the  debond  region  with  length  /j,  there  exists  a  frictional  shear 

stress  X  which  is  proportional  to  the  normal  pressure  p  at  the  interface. 
Through  the  detail  derivation  in  Appendix  A,  the  stress  fields  in  the 
debond  region  are  evaluated  as  follows: 


=  c[exp(y9x^  /r)  - 1]  +  CT^ 

(11) 

(12) 

T  =  -(/3c/2)exp()3.r^/r) 

(13) 

p  =  {Pc/2iu)exp[pxJr) 

(14) 

where 

c=  T}-'. 

n  and  ^  are  expressed  in  terms  of  Young’s  modulus  E^,  E„,  Poisson  ratio  v^, 
v„  and  fiber  volume  fraction  /  in  Appendix  A.  p  is  the  friction  coefficient. 
Afl  is  the  residual  thermal  strain  defined  as 

=  (15) 

where  and  a„  are  the  thermal  expansion  coefficients  of  fiber  and  matrix 
respectively,  A7(<0)  is  the  temperature  change  from  the  stress-free 
temperature  of  matrix  to  the  operating  temperature.  The  interface  is 
subject  to  compression  only  when  p  is  negative.  According  to  Eq.(14), 
following  condition  must  be  satisfied 
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4i3+(v^(T^y'£^)<0. 

Through  the  whole  paper,  condition  (16)  is  always  assumed.  Now,  one  can 
calculate  the  crack  opening  at  the  bridged  domain  due  to  the  elastic 
extension  of  the  fiber  out  of  the  matrix.  The  crack  opening  displacement  is 
defined  as 


One  can  integrate  and  e„  given  in  Appendix  B  from  0  to  such  that 
and  u„  at  the  crack  plane  are  obtained.  The  crack  opening  displacement  is 
found  to  be 


(18) 


From  the  eigenstrain  approach,  the  crack  opening  is  defined  as 

5  =  2c„ae  =4rfi^Y-^'|(T  (19) 

at  the  bridged  domain  0^.  By  setting  Eq.(18)  equal  to  Eq.(19),  a  is 
determined  and  expressed  as 


(20) 

It  is  clear  that  a  =  0  as  4=0  (no  sliding),  a  =  1  as  4  (complete  sliding 
for  a  continuous  fiber). 

2.3  Analysis  of  debond  length 

The  debond  length  a  be  determined  by  using  the  energy  balance 
approe  ,h.  If  the  debond  crack  extends  by  an  incremental  amount  dl  at 
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constant  o^,  several  different  kinds  of  energy  change  will  occur  and  reach 
to  the  equilibrium  state.  The  strain  energy  is  increased  by  dE,,  work  is 
done  by  traction  stress  with  magnitude  dw^,  frictional  energy  dw^  is 
dissipated,  debonding  energy  dw^  is  absorbed  during  interfacial  debonding. 
The  energy  balance  condition,  at  /  =  /j,  can  be  formulated  as  below: 

dw^-dE  -dw.-dw.-Q.  (21) 

T  s  d  f 

Dividing  each  term  by  the  incremental  debond  length  dl,  one  can  change 
Eq.(21)  into  the  following  expression: 
dw  dE  dw^  dw 

— L - 1  =  — d_  +  — L.  (22) 

dl  dl  dl  dl 

Each  energy  term  is  derived  in  Appendix  B.  For  small  exponent  index 
p/j/r,  Eq.(22)  is  used  to  obtain  the  approximate  value  of  I/. 
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Here,  p  is  defined  in  Appendix  B.  Furthermore,  substituting  Eq.(23)  into 
Eq.(20),  one  can  obtain  the  approximate  value  of  a  for  small  p/^/r  as 


a  =  l-<l  + 


/.(.j  E 


L-llfic-ACi 


In  Eq.(24),  a  is  a  function  of  bridging  stress  Oj  and  several  material 


constants. 


2.4  Energy  calculation  for  a  bridged  crack 


ixterni 


In  order  to  calculate  the  energy  release  rate  during  the  crack 
growth,  the  Gibbs  free  energy  of  the  composite  under  an  applied  stress  is 
evaluated  when  the  configuration  of  Fig.  1  is  simulated  by  inclusion  Qq  and 

inclusions  Qx.  Also,  the  matrix  is  approximated  by  the  homogeneous 
medium  with  the  composite  moduli  jl  and  v.  According  to  the  definition 

(Mura,  1987),  the  Gibbs  free  energy  is 

"  "I  k  '  k 

where  is  given  by  Eq.(5)  and  Eq.(6),  ejj  is  given  by  Eq.(3)  and  Eq.(4). 
Substituting  them  into  Eq.(25),  F  is  evaluated  as 

4(l-v)a’(<T  f  l+(l-a^)-/ 
f  =  — AI - - Lr — .  (26) 
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By  definition,  the  energy  release  rate  of  a  crack  is  written  as 

G  =  — ^ 
dna^ 

which  is  evaluated  below: 
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where 

C„=2(l- v)(er^)^<j/;r^. 

Go  is  the  energy  release  rate  of  a  penny-shaped  crack  in  the  homogeneous 
medium  without  fibers.  (3a/3a)  can  be  evaluated  by  the  following 
procedures:  First,  one  can  substitute  Eq.(24)  into  Eq.(lO)  and  obtain 
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(29) 

Second,  {dafda)  is  equal  into  {da/ do j)/ {da/ da j).  Using  Eq.(24)  and  Eq.(29)  to 
obtain  {da/dOf)  and  {da/dOj)  individually,  one  can  calculate  {da/da) 
analytically  or  numerically. 


2.4.2  Dissipation  energy  and  debonding  energy 

The  crack  opening  in  Q,  is  achieved  when  the  fiber  extends  out  of 

the  matrix.  This  is  an  inelastic  deformation,  confirming  to  the  definition  of 
eigenstrain  ae^.  Before  the  fiber  begins  to  slide,  interfacial  debonding 

must  initiate  at  first.  Therefore,  dissipation  energy  due  to  the  work 
against  friction  and  debonding  energy  needed  to  create  new  surface  at  the 
interface  can  contribute  to  the  increase  of  crack  resistance.  One  can 
define  the  total  dissipation  energy  and  the  total  debonding  energy  as 
follows: 


(30) 

(31) 


where 


8 


r^' 

N  is  the  amount  of  fibers  bridging  the  matrix  crack.  2N  represents  the 
amount  of  debond  lengths  on  the  both  sides  of  the  matrix  crack.  Using 
Eq.(B7)and  Eq.(B8),  one  obtains 


(32) 

(33) 


where 


/(ctJ  =  |pc^  -  +  cE^AQ^.  (34) 

One  can  do  {dajdcSj)  and  (^/(aj/^cr^.)  using  Eq.(23),  (29)  and  (34). 

Then,  {dl^/da)  and  (9/(a  )/3a)  are  calculated  by  applying  the  chain  rules  of 
differentiation.  The  total  crack  resistance  is  defined  as 


dW  dW 

(35) 

dna^  dna^ 

where  is  the  inherent  fracture  toughness  of  the  matrix,  the  second  and 
third  terms  on  the  right-hand  side  of  Eq.(35)  are  the  additional  crack 
resistances  due  to  frictional  sliding  and  interfacial  debonding, 
respectively. 


3.  DISCUSSION 

In  this  section,  we  are  going  to  discuss  how  the  thermal  residual 
effect  influences  fracture  mechanisms  such  as  interfacial  debonding,  fiber 
failure  and  matrix  cracking.  Following  parameters  are  used  for  numerical 
calculation:  E^  =  5S0GPa,  £„  =  410 GPa,  v^=0.23,  v„=0.19,  a^  =  4.8x10^, 

a„  =  8.9  X 10^,  /  =  0.2 ,  p  =  0.1,  r  =  0.4  pm,  C,  =  30  Jlm^ . 

3.1  Fracture  mechanisms 

3.1.1  Interfacial  debonding 


9 


It  follows  from  Eq.(23)  that  /j->0  as  approaches  to  I. 
Accordingly,  Z  is  a  threshold  stress  for  debonding  initiation.  One  can 
obtain  E  by  solving  the  following  equation: 

2E  G 

g((Tr) - ^  =  0  (36) 


which  has  the  solution 

1 


(7=1  = 

^  2D, 


+  D, 


SE^G 


(37) 


In  Fig. 3,  it  clearly  shows  that  E  increases  as  G,  increases  for  the  same  AT. 
Also,  for  the  same  G,,  E  decreases  as  \AT\  increases.  The  reason  is  before 
debonding,  when  \AT\  becomes  larger,  the  fiber  is  under  higher  residual 
compression  in  direction  and  the  surrounding  matrix  is  under  higher 
residual  tension.  After  debonding,  the  fiber  tends  to  extend  and  matrix 
tends  to  shrink  naturally  because  of  relaxation  of  axial  residual  stresses. 
So,  debonding  initiation  is  easier  to  occur.  The  same  tendency  was  also 
predicted  by  Hsueh  (1991)  for  the  case  of  isotropic  thermal  expansion 
coefficients  of  the  fiber.  Eq.(37)  reconfirms  the  previous  conclusion  that  E 
is  lower  as  r  is  larger  (Yang  et  al.,  1991).  Replacing  cr,-  in  Eq.(lO)  by  E  and 
letting  a  =  0  for  0,  one  obtains 


rK/^-/r 


(38) 


where  Z  is  a  function  of  AT  and  G^.  For  given  AT  and  G,  ,  curve  of  Eq.(38), 
debonding  initiation,  is  shown  in  Fig.4.  As  ajr  » 1 ,  has  an  asymptotic 
value: 

(39) 

If  the  Poisson  effect  and  thermal  expansion  mismatch  are  ignored  by 
setting  v^=0  and  AT  =  0,  one  obtains  and  D,  =0.  The  threshold 

stress  is  expressed  as 


E  = 


\4EG 

f  S 


(40) 


which  is  exactly  identical  to  that  of  the  frictionless  case  studied  by 
Outwater  and  Murphy  (1969). 


3.1.2  Fibgr  feilvrg 

It  is  assumed  that  fibers  have  the  averaged  uniform  strength  5. 
Fiber  failure  occurs  when  bridging  stress  reaches  fiber  strength,  i.e. 

(41) 


Substituting  Eq.(41)  into  Eq.(lO),  one  obtains 


(42) 


where 


a 

s 


'  0 
'  a(S) 


I<S< -E^AQ jv ^ ' 


Replacing  cSj  by  5,  a(s)  is  obtained  from  Eq.(24).  Fig.5  is  the  plot  of  Eq.(42) 
for  two  different  AT  and  5  .  Obviously,  as  S  increases,  must  increase  to 
cause  fiber  failure  for  the  same  crack  length  at  the  same  AT .  It  is  worth 
noting  that  there  is  a  critical  value  of  strength  to  separate  the  opposite 
effects  of  AT  through  in  Eq.(42).  As  S<S^,  larger  \AT\  (larger  a^)  causes 
fiber  failure  at  larger  applied  stress.  As  S  >  S^,  larger  \AT\  (  smaller  a,) 
results  in  fiber  failure  at  smaller  applied  stress.  Only  when  S  =  S^,  these 

two  AT  has  the  same  effect  on  fiber  failure  mechanism.  Fig.6  shows  the 
relation  between  a,  and  S  for  two  different  AT.  From  Eq.(42),  one  obtains 
the  asymptotic  value  of  as  a/r  » 1 : 

0-^=yS.  (43) 


3.1.3  Matrix  cracking 

During  crack  growth  equilibrium  between  energy  release  rate 
(crack  driving  force)  and  the  crack  resistance  can  be  expressed  by 

G  =  R  (44) 

Substituting  Eq.(28)  and  Eq.(35)  into  Eq.(44),  one  obtains  the  critical  value 
of  (Jj  causing  matrix  cracking.  Oj  is  expressed  as  a  function  of  and  a. 
This  function  can  be  substituted  into  Eq.(29)  to  plot  the  curve  of  matrix 
cracking  in  Fig.7.  One  can  also  prove  numerically  that 

dG  dR 

- - - r 

dna  dna 


1  1 


(45) 


which  is  satisfied  for  any  given  critical  applied  stress  Therefore,  the 
matrix  crack  always  grows  unstably  at  As  a/r»l,  the  asymptotic 
values  of  G  and  R  are  written  below: 


4(l-w);i 

(^) 

( l  +  a^ 

37tfif  ' 

U-aJ 

/?  =G  + 

oo  c 


^  dW  ^ 

dis. 


dna" 


+ 


deh. 


J 


dna" 


J 


(46) 

(47) 


where 


dis. 


dna" 


J 


■f'hfc/')’ 


/ 


Ow  ^  ,  , 


dna‘ 


J 


(48) 

(49) 


According  to  Eq.(29),  as  afr  <»,  approaches  to  tr^//.  Therefore,  one  can 
replace  Oj  by  a^/f  to  obtain  the  asymptotic  expressions  of  4,  a  and  /(ct,) 
from  Eq,(23),  (24)  and  (34),  respectively.  Generally,  fracture  criterion, 
yields  an  algebraic  equation  for  critical  applied  stress  cr“.  For 
complete  bridging  case  (a  =  0),  G  and  R  can  be  evaluated  analytically  as 
follows: 


G  =G  3-(-2/a/r 

“  =  '  ^3{\  +  fa/rf 


R 


a  =  0 


The  as)nnptotic  values  of  G„^o  ^a=o»  respectively,  are  given  by 


(50) 

(51) 

(52) 

(53) 


Letting  Eq.(52)  be  equal  to  Eq.(53),  the  asymptotic  value  of  critical  applied 
stress  is  given  as 


1  2 


(54) 


For  the  given  material  parameters  in  the  section  of  discussion,  Fig. 7 
shows  the  curves  of  matrix  cracking  for  two  different  AT .  For  the  same 
crack  length,  smaller  \AT\  makes  matrix  cracking  occur  at  slightly  lower 

applied  stress.  One  can  not  see  much  difference  between  these  two  curves. 
However,  if  one  uses  the  curves  of  Eq.(38),  Eq.(42)  and  matrix  cracking  to 
construct  the  fracture  map  as  a  (a/r  V5  a^)  diagram,  AT  will  show  its 
influence  on  the  crack  extension  paths  and  final  failure  mode  of  a 
composite.  Fig.8  and  Fig.9  are  the  fracture  maps  for  larger  and  smaller 
\AT\,  respectively.  For  a  large  crack,  path  (A)  shown  in  Fig.8  hits  the  curve 

of  interfacial  debonding  first.  Further  increase  of  to  point  A  will  result 
in  matrix  cracking  through  whole  matrix  with  no  fiber  failure.  It  is  a  non- 
catastrophic  type  of  failure  and  a  precursor  of  multiple  matrix  cracking 
causing  larger  nonlinear  strain.  The  final  composite  strength  can  be 
defined  as  fS  from  Eq.(43).  For  smaller  \AT\  (smaller  residual  effect),  path 

(A)  shown  in  Fig.9  also  passes  the  curve  of  interfacial  debonding  but 
increasing  causes  fiber  failure  directly  without  matrix  cracking.  It  is  a 
catastrophic  type  of  failure.  The  composite  strength  is  determined  by 
point  A  obtained  from  Eq.(42).  For  intermediate  and  small  cracks,  path  (B) 
and  (C)  in  Fig.8  indicate  that  debonding,  matrix  cracking  and  fiber  failure 
occur  in  sequence.  The  composite  fails  catastrophically  at  stress  level  of 
point  B  and  C  for  intermediate  and  small  cracks,  respectively.  Composite 
strength  is  predicted  by  substituting  initial  crack  length  into  Eq.(44).  In 
Fig.9,  there  is  a  critical  point  O  associated  with  initial  crack  length  Oq.  For 
an  intermediate  crack,  a>aQ,  path  (B)  passes  through  the  curve  of 
debonding  initiation  and  reaches  the  curve  of  fiber  failure.  It  implies  that 
the  whole  composite  fails  catastrophically  at  point  B  which  is  determined 
by  Eq.(42).  For  a  small  crack,  a<ao,  path  (C)  is  similar  to  that  in  Fig.8. 

From  the  comparison  between  Fig.8  and  Fig.9,  one  can  discuss  how  AT 
changes  the  fracture  mechanisms  of  a  composite  as  well  as  the  final 
composite  strength  for  a  given  initial  crack  length.  Generally,  the 
fracture  map  is  also  useful  to  discuss  the  effects  of  other  parameters  such 
as  G^,  fi,  f  and  r  on  the  fracture  mechanisms  of  a  composite. 


3.2  Fracture  transition  for  a  large  crack  (a/r»l) 


Two  categories  of  fracture  transitions  will  be  discussed  in  this 
section.  The  critical  conditions  of  transitions  are  also  defined. 

3.2.1  Complete  bridering  versus  interfacial  debonding 

Using  Eq.(39)  and  Eq.(54),  one  can  define  a  characteristic  constant: 

4(1- vV/X^ 

K  =  -^ - ^ - .  (55) 

As  fc  <  1 ,  debonding  occurs  at  the  interface  before  the  matrix  crack  extends. 
As  ol,  debonding  can  not  occur  during  the  matrix  cracking  through 
whole  matrix.  Therefore,  k  =  1,  it  is  the  critical  condition  for  interfacial 
debonding. 

3^.2  Catastrophic  failure  versus  non-catastrophic  failure 

First,  K<1,  non-catastrophic  failure  is  predicted  under  the 
condition: 

fl<cr  <  fS.  (56) 

Using  inequality  (56),  one  can  find  the  range  of  certain  parameter  such  as 
47,  /,  r  et  al.  to  facilitate  non-catastrophic  type  of  failure.  Second,  )f>l, 
the  following  inequality 

<fS  (57) 

ought  to  be  satisfied  for  non-catastrophic  failure  of  a  composite.  Inecuality 
(57)  can  be  rewritten  as 

<  fS.  (58) 

It  is  noticeable  that  AT  does  not  involve  in  inequality  (58).  It  means  if 
interfacial  debonding  does  not  occur,  residual  thermal  effect  has  no 
influence  on  the  transition  between  catastrophic  and  non-catastrophic 
failure.  Take  the  crack  paths  (A)  in  Fig.8  and  Fig.9  as  examples.  For  Fig.8, 
x’  =  0.28.  For  Fig.9,  k-  =  0.50.  Beth  values  of  k  are  smaller  than  one  so  that 
interfacial  debonding  always  initiates  before  the  matrix  crack  extends. 
For  Fig.8,  /Z  =  0.82GPa,  a"=1.06GPa,  and  /S  =  1.10 GPa  satisfy  inequality 
(56).  For  Fig.9,  /r=1.08GPa,  <t“  =  l.C3GPa  and  fS=\.\QGPa  do  not  satisfy 


4(1-  v)r 


inequality  (56).  Therefore,  non-catastrophic  failure  is  only  expected  for 
path  (A)  in  Fig.8. 


4.  CONCLUSION 

The  thermal  residual  effect  on  the  fracture  mechanisms  including 
interfacial  debonding,  matrix  cracking  and  fiber  failure  are  investigated 
in  this  paper.  Combining  Eshelby’s  inclusion  method  and  fracture  criteria 
of  fracture  mechanics,  the  fracture  mechanisms  are  governed  by  Eq.(38) 
for  debonding,  Eq.(42)  for  fiber  failure  and  Eq.(44)  for  matrix  cracking. 
The  fracture  maps  like  Fig.8  and  Fig.9  are  used  to  analyze  the  crack 
extension  paths  and  final  composite  strength.  Comparing  Fig.8  with  Fig.9, 
one  can  conclude  that  AT  has  the  strong  influence  on  the  crack  extension 
paths  and  final  composite  strength.  For  the  given  material  parameters  in 
the  section  of  discussion,  decreasing  \AT\  (decreasing  thermal  residual 

effect)  undergoes  the  transition  from  non-catastrophic  type  to  catastrophic 
type  of  failure  for  a  large  crack  in  a  fiber-reinforced  brittle  matrix 
composite.  In  general,  the  conditions  of  fracture  transitions  for  a  large 
crack  are  defined  in  section  3.2.  These  conditions  may  provide  some  useful 
guidelines  to  prevent  undesirable  fracture  mechanisms  when  a  composite 
is  designed. 
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Appendix  A 


Consider  the  schematic  of  Fig. 2.  Assuming  zero  radial  stress  on  the 
outside  of  the  composite  cylinder  and  using  a  modified  shear  lag  model 
(Sigl  and  Evans,  1989),  the  stress  fields  satisfy  the  Lame’  relationships: 

<^i=<^L=p  •  p<o 
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where  is  the  radius  of  the  cylinder  composite,  r  is  the  fiber  radius,  r  is 

the  distance  from  the  center  of  the  fiber  to  the  point  in  the  matrix  and 
superscripts  /  and  m  denote  fiber  and  matrix,  respectively.  The 
corresponding  strains  including  thermal  residual  strains  are 
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The  requirement  of  continuity  of  radial  displacements  at  the  interface 
r  =  r: 


pf  _  pm 
66  66 

jdelds  the  relation  between  normal  pressure  and  axial  stresses  as 
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where 
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+  1- V 


Assume  that  the  interfacial  friction  is  governed  by  Coulomb’s  friction  law 
T  =  -fip  (All) 

where  p.  is  a  constant  friction  coefficient.  Moreover,  the  equilibrium 
equations  in  the  direction  must  be  satisfied  within  the  debond  region, 
0<X}<1^.  They  are  written  as 


dcL 
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Combining  Eqs.  (AlO),  (All),  (A12)  and  (A13),  one  obtains  the  following 
differential  equation 

da{ 


where 


I3  =  2nii(  ,  »7=  w  ±. 
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When  boundary  conHitions 

<T3"3  =  0  ar^3=0 

are  applied  into  Eq.(A14),  the  solutions  are  shown  as  below: 
<^3^3  =  4exp(/3x^ /r)  -  ll  +  (j^ 


r-=-c[exp(PxJr)-\]\^-^ 
■  =  -{Pc/2)exp(PxJr'j 
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where 


Appendix  B 


Using  the  stresses  and  strains  obtained  in  Appendix  A,  one  can 
calculate  strain  energy  energy  dissipation  due  to  frictional  sliding, 

debonding  energy  absorbed  by  the  interface  and  work  Wj  done  by  the 
traction  stress  The  elastic  strain  energies  in  the  fiber  and  matrix  are 
associated  with  and  £33.  This  estimation  gives 
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The  debonding  energy  is 

w.=2KrlG  (B2) 

a  s 

where  is  the  fracture  toughness  of  the  interface. 

The  energy  dissipation  due  to  frictional  sliding  is 

w^=2jTrl^T{u^-u^)dx^.  (B3) 

The  work  done  by  the  traction  stress  is 

=  7tr^aJu-u  I  (B4) 

^  ^1x^=0 

where  and  u„  are  obtained  by  integrating  £«  and  in  Appendix  A, 
respectively.  For  small  {Pl/r),  each  energy  term  can  be  simplified  as 
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where  are  defined  in  Appendix  A. 


FIGURE  LEGENDS 


Fig.l  Inclusion  model  for  the  crack  and  the  cross  sections  of 

fibers  at  the  crack. 

Fig. 2  The  rod-cylinder  model:  bridging  stress  is  balanced  by 
the  frictional  stress  x  caused  by  the  normal  pressure  P 
acting  at  the  debond  interface. 

Fig. 3  Threshold  stress  Z  for  debonding  plotted  as  a  function  of 
interfacial  fracture  toughness  for  two  different 

temperature  change  AT. 

Fig. 4  Curves  of  debonding  initiation  plotted  in  a  crack  length  vs 

applied  stress  diagram. 

Fig. 5  Curves  of  fiber  failure  plotted  in  a  crack  length  vs  applied 
stress  diagram. 

Fig. 6  Sliding  parameter  a,  plotted  as  a  function  of  fiber  strength. 

Fig. 7  Curves  of  matrix  cracking  plotted  in  a  crack  length  vs 
applied  stress  diagram. 

Fig. 8  Fracture  map  for  zl7'  =  -980V.  Dash  lines  (A),  (B)  and  (C), 

respectively,  represent  the  extension  paths  for  large, 
intermediate  and  small  cracks. 

Fig. 9  Fracture  map  for  zlT’  =  -490V.  Critical  point  O  is  associated 

with  initial  crack  length  a^. 
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Abstract 

In  single  fiber  pullout  experiments  typical  configurations  of  loading  include  support 
of  the  matrix  on  the  side  of  the  protruding  fiber  and  support  of  the  matrix  firom  the 
opposite  end.  Single  fiber  pullout  experiments  have  been  conducted  using  these  two 
configurations,  and  the  experimental  results  are  compared.  The  comparison  shows 
that  the  difference  in  the  measured  fiber  pullout  fidctional  stresses  between  the  two 
loading  configurations  is  minimal  for  small  embedded  lengths.  However,  the 
difference  increases  with  the  fiber  embedded  length.  Theoretical  analysis  predicts  the 
same  trends  with  loading  configuration  change. 

1.  Introduction 

Fiber-reinforced  composites  offer  the  potential  to  meet  the  demands  of  materials 
that  are  light,  strong  and  reliable.  The  mechanical  properties  of  these  materials 
are  sensitive  to  the  fiber/matrix  interfadal  properties,  in  addition  to  the 
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intrinsic  properties  of  the  fibers  and  the  matrices  [1,2].  Recognizing  the  importance 
of  interfadal  properties  has  resxilted  in  numerous  studies  to  characterize  the 
debonding  and  fidctional  sliding  phenomena  [3-8].  Fiber  pullout  and  push-down 
(indentation)  are  the  two  most  firequently  used  experimental  techniques  for 
measuring  the  interfacial  debonding  stress  and  firictional  stress  in  fiber  composites. 
Detailed  analyses  of  these  two  different  methods  of  loading  of  the  fiber  show 
substantial  differences  between  fiber  pullout  and  push-down  due  to  the  Poisson  effect 
in  the  fiber  [9].  In  fiber  pullout  testing,  different  methods  of  matrix  loading  have  also 
been  used,  including  support  of  the  matrix  opposite  the  side  of  the  protruding  fiber 
and  support  of  the  matrix  on  the  side  of  the  pi:otruding  fiber  ^ 
as  illustrated  in  Fig.  1.  Such  a  difference  in  loading  methods  for 
fiber  pullout  may  also  restdt  in  different  measured  interfadal  properties  due  to 
different  boxmdary  conditions.  The  objective  of  this  woric  is  to  compare  the  frictional 
pullout  stresses  obtained  using  these  two  methods  of  matrix  loading  during  pullout. 

2.  Experimental  procedure 

The  specimen  geometry  and  loading  fixtures  are  shown  in  Fig.  2  for  loading 
configurations  A  and  B,  respectively.  In  both  fixtures,  a  self-aligning  mechanism  was 
used  to  ensure  good  aUgnment  for  fiber  pullout.  In  the  specimen  geometry  for  loading 
configuration  A,  the  fiber  was  pre-notched  so  that  the  fiber  failed  at  the  notch  upon 
loading  and  was  subsequently  pulled  out. 

an 

The  single  fiber  composite  consisted  of  steel  piano  wire  acrylic^ 

matrix.  The  acrylic  fails  in  a  brittle  manner  at  room  temperature.  The  steel  wire 
has  a  diameter  of  0.23  mm,  and  has  a  yield  strength  of  3  GPa,  and  fails  with  a  very 
small  degree  of  yielding.  Debonding  and  pullout  of  the  steel  wire  took  place  at  stresses 


5  Acrylite  FF,  GYRO  Industries,  Mt.  Arlington,  NJ. 
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well  below  the  yield  point.  The  steel  wire  was  cleaned  with  acetone  and  methyl 
alcohol,  and  the  acrylic  was  cleaned  with  methyl  alcohol  before  being  used  for 
specimen  preparation.  The  pullout  specimens  were  fabricated  by  sandwiching  the 
steel  wire  between  two  acrylic  plates  and  hot  pressing  at  155±2°C.  The  pullout 
specimens  then  were  machined  from  the  hot  pressed  samples. 

Pullout  experiments  were  performed  itaing  an  Instron  universal  testing  machine^ 
at  a  constant  crosshead  speed  of  0.5  mm/min.  The  load  and  crosshead  displacement 
were  recorded  using  a  computer  data  acquisition  system.  Stick-slip  behavior  was 
observed  as  the  fiber  was  pulled  out  the  matrix.  The  pullout  load  was  taken  from  the 
peak  load  of  each  stick-slip  section.  The  corresponding  fiber  embedded  length  was 
determined  from  the  load-displacement  curve  after  subtracting  the  displacement  due 
to  the  load  train  compliance. 

4.  Results  and  discussion 

The  results  of  frictional  pullout  stress,  <^,  vs.  fiber  embedded  length,  ft)!*  the  two 
loading  configurations  are  shown  in  Fig.  3  for  fiber  embedded  lengths  ranging  firom  2  to 
12  mm.  Results  are  shown  for  five  pxillout  tests  for  each  configuration.  The  dashed 
line  and  the  solid  line  are  the  curve  fits  to  the  data  using  fimction  (^  =  Co  [1  -  exp(ci  /«)] 
for  results  of  loading  configurations  A  and  B  where  Co  and  Ci  are  constants.  While 
there  is  a  scatter  in  the  data,  a  comparison  of  the  curves  shows  a  difference  in  pullout 
stresses  in  the  upper  range  of  the  fiber  embedded  lengths  examined,  with  the  pullout 
stress  for  loading  configuration  B  lower  than  that  for  loading  configuration  A.  The 
difference  in  the  pullout  stresses  between  the  two  loading  configurations  increases 
with  the  fiber  embedded  length.  It  is  expected  that  the  difference  would  further 
increase  with  the  fiber  embedded  length. 


^  Instron  Model  1125,  Instron,  Corp.,  Canton,  MA. 
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The  difference  in  the  pullout  stress  between  the  two  loading  configurations  can  be 
understood  physically  as  follows.  In  loading  configuration  B,  the  matrix  is  in  axial 
compression.  Due  to  the  Poisson  effect,  axial  compression  increases  the  diameter  of 
the  matrix  hole  that  holds  the  fiber.  Therefore,  a  lower  puUout  stress  is  expected.  In 
loading  configuration  A,  axial  tension  in  the  matrix  decreases  the  diameter  of  the 
matrix  hole  that  holds  the  fiber,  and  thus,  a  greater  pullout  stress  is  anticipated.  The 
Poisson  effect  in  the  matrix  as  observed  is  considerably  smaller  than  that  in  the  fiber 
because  the  axial  stress  in  the  matrix  is  significantly  lower  than  that  in  the  fiber  as 
the  matrix  size  is  many  times  that  of  the  fiber  in  typical  single  fiber  pullout  tests. 

The  experimental  results  of  the  present  study  are  compared  with  the  the 
theoretical  analysis  by  Hsueh  [10]  (see  Appendix  A  for  more  details  of  the  analysis). 
It  is  assumed  that  both  the  matrix  and  the  fiber  display  linear  elastic  behavior,  and 
that  the  interfadal  pullout  may  be  described  by  Coulomb  fiiction.  In  order  to 
compare  experiment  with  theory ,  two  parameters  which  describe  the  Coulomb 
firiction  at  the  interface  are  needed:  the  normal  residual  stress  (cTc)  cuid  a  coefficient  of 
fiiction  (p).  This  two  parameters  were  obtained  experimentaUy  using  a  technique 

described  by  Takaku  and  Arridge  [11].  In  this  experiment,  the  dependence  of  the 
pullout  stress  (c^)  on  the  initial  tensile  stress  on  the  fiber  was  determined  and 

analyzed  using  the  following  equation  to  obtained  residual  normal  compressive  stress 
and  the  coefficient  of  fiiction: 

Op  =  ( ^  -  g/ )  [l  -  exp  ( -  —  ^  )]  (1) 

where  le  is  the  fiber  embedded  length,  a  is  the  radius  of  the  fiber, 
k  =  (J&mVfVDEf  (l+Vm)],  Ejn,  E{,  Vm,  and  vrare  the  Yoimg's  moduli  and  the  Poisson's 
ratios  of  the  matrix  and  fiber.  For  the  steel  wire-acrylic  system,  the  residual  normal 
compressive  stress  was  determined  to  be  8  MPa  with  a  coefficient  of  friction  p  =  0.5. 
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With  the  measiired  normal  compressive  stress  and  the  coefficient  of  friction,  the 
predicted  pullout  stresses  for  the  two  loading  configurations  are  computed  as  a 
function  of  the  embedded  length,  and  are  compared  with  experimental  results  in 
Fig.  4(a).  Although  the  experimental  and  theoretical  results  do  not  agree  well,  the 
effect  of  the  loading  configurations  demonstrates  the  same  trend,  that  is,  the  pullout 
stress  determined  via  loading  configuration  A  is  greater  than  that  fi*om  loading 
configuration  B  for  a  given  embedded  length. 

An  alternative  way  to  compare  the  experimental  and  theoretical  results  is  to  fit 
the  experimental  data  for  one  loading  configuration  to  obtain  the  normal  compressive 
stress  and  the  coefficient  of  firiction  and  to  compute  the  expected  behavior  for  the 
second  loading  configuration.  These  results  may  then  be  compared  with  the 
experimental  results  for  the  latter  loading  configuration.  Such  a  comparison  is  shown 
in  Fig.  4(b),  in  which  the  experimental  results  for  loading  configuration  A  is  fit  into 
eqn.  (1)  in  Ref  10  with  the  coefficients  Ai  and  A2  in  ^pendix  A,  giving  a  normal 
compressive  stress  of  2.2  MPa  and  a  fnction  coefficient  of  2.3.  The  theoretical 
prediction  for  loading  configuration  B  is  computed  using  eqn.  (8)  in  Ref.  10  and  is 
plotted  along  with  the  experimental  results  for  the  same  loading  configuration.  The 
comparison  now  shows  good  qualitative  agreement  between  the  experiment  and  the 
theory.  A  difficulty  with  this  approach  is  that  the  coefficient  of  friction  obtained  is 
unreasonably  high.  Furthermore,  the  obtained  values  of  the  coefficient  of  firiction  and 
normal  compressive  stress  are  inconsistent  with  those  obtained  using  Takaku  and 
Arridge  technique. 

In  order  to  gain  a  better  understanding  of  the  problem,  the  fiber  surface  and  the 
matrix  surface  at  the  interface  before  and  after  fiber  sliding  were  examined  using 
scanning  electron  microscopy  [12].  The  observation  shows  that  the  matrix  siuface 
at  the  interface  changes  significantly  after  fiber  sliding.  Loose 

debris  was  observed  on  the  interface  matrix  surface  after  sliding.  At  the  same  time, 
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the  matrix  surface  was  roughened.  Therefore,  the  coefficient  of  firiction  at  the 
interface  changes  with  sliding  and  with  sliding  distance  of  fiber  pullout.  The  normal 
compressive  stress  might  have  changed  also  due  to  the  matrix  wearing.  These 
changes  in  the  coefficient  of  firiction  and  normal  compressive  stress  with  aliHing 
distance  could  explain  the  poor  agreement  between  the  experiment  and  theory,  the 
latter  of  which  assumes  a  constant  coefficient  of  friction  and  normal  compressive 
stress.  This  observation  has  a  strong  implication  for  detracting  the  interfadal 
properties  firom  single  fiber  pullout  tests:  misleading  values  can  result  firom  fitting  the 
data  to  models  which  assume  constant  coefficient  of  friction  and  normal  compressive 
stress  if  in  reality,  these  parameters  change  with  fiber  puUout. 

5.  Summary 

The  fiber  fiictional  pullout  stresses  of  two  matrix  loading  configurations  have  been 
compared  experimentally.  The  difference  of  results  between  the  two  matrix  loading 
configuration  is  minimal  for  short  fiber  embedded  lengths.  The  difference  increases 
with  fiber  embedded  length,  indicating  the  need  to  consider  the  matrix  loading 
condition  for  fiber  pullout  with  long  fiber  embedded  lengths.  A  comparison  between 
theory  and  experiment  was  complicated  by  nhanging  coefficient  of  fidetion  and 
probably  changing  normal  compressive  stress  with  fiber  pullout.  However,  the 
comparison  shows  a  qualitative  agreement,  that  is,  the  pullout  stress  for  loading  the 
matrix  face  opposite  the  protruding  fiber  is  greater  than  that  obtained  firom 
supporting  the  matrix  on  the  side  of  the  protruding  fiber. 
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Appendix  A.  Theoretical  analysia  of  fiber  Motional  pullout 

Hsueh  [10]  analyzed  the  two  loading  configurations  (Fig.  1)  for  pullout  tests  using 
an  idealized  model  for  a  single  fiber  pullout  consisting  of  a  fiber  located  in  the  center  of 
a  co-axial  cylindrical  matrix.  The  pullout  stresses  for  the  two  loading  configurations 
are  given  by  Hsueh  in  Ref.  10.  Due  to  the  length  of  the  equations,  they  are  not 
reproduced  here. 

In  Hsueh's  work  [10],  an  approximation  of  the  shear  stress  distribution  in  the 
matrix  results  in  a  finite  value  of  shear  stress  at  the  outer  surface  of  the  matrix. 
Hsueh  [13]  improved  the  original  analysis  by  introducing  an  alternative  shear  stress 
distribution  in  the  matrix  which  satisfies  the  boundary  conditions  of  zero  tractions  on 
the  cylindrical  surface  of  the  matrix.  The  results  of  the  pullout  stress  are  the  same 
as  those  given  in  the  original  work  [10],  with  the  exception  of  changes  in  coefficient  Ai 
and  A2,  which  are  given  as  [13]: 


Al  =! 


a(l-6*/o*)D 


2  Vn,  (1  +  vj  L  ^  ln(6/a)- - - J 


(Al) 


Ao  = 


(l-6"/a")(g„Vf/gfV„)-l 


o  M  .  6^  r.o .  (6^-Q=^)(3&^-a")i 

2  /i  Vm  ( 1  +  Vn,)  ^2 ^2  L  ^  ln(i>/fl)  —  ^  j^2  J 


(A2) 


46* 


The  effects  of  this  change  on  the  pullout  stress  is  very  small.  However,  the  new 
solutions  are  used  because  they  satisfy  the  boundary  conditions  of  free  tractions  on 
the  outer  cylindrical  surface  of  the  matrix. 


9 


Figure  Captions 

Fig-  i-  Schematic  of  the  two  matrix  loading  configurations  used  in  fiber  pullout  tests: 
(a)  support  of  the  matrix  opposite  the  side  of  the  protruding  fiber  and  (b) 
support  of  the  matrix  on  the  side  of  the  protruding  fiber. 

Fig.  2.  Schematic  of  specimen  geometry  and  loading  fixture  for:  (a)  loading 
configuration  A  and  (b)  loading  configuration  B. 

Fig.  3.  Experimental  results  comparing  fiictional  pullout  stress  for  the  two  different 
loading  configurations  examined. 

Fig.  4.  (a)  A  comparison  of  experimental  results  and  theoretical  predictions  based  on 
an  independent  measurement  of  the  normal  compressive  stress  and 
coefficient  of  friction;  (b)  A  comparison  of  experimental  results  and 
theoretical  prediction  for  configuration  B.  The  normal  compressive  stress 
and  coefficient  of  firiction  were  obtained  by  curve  fitting  the  experimental 
results  to  Hsueh's  theory  [10,13]  for  configuration  A. 


10 


Fig.  1 


o 


screw 
connection 

/ 


IIIB 


Embedded  Length,  /,  (mm) 


IG186(LO:GO) 


Pig.  4 


IG186(L1:G1,G3) 


APPENDIX  5 


Inverse  Problems  in  Linear  Elasticity 

by 


Tatsuhito  Koya 
and 

Toshio  Mura 


Department  of  Civil  Engineering 
Northwestern  University 
Evanston,  Illinois  60201 


Abstract 

Some  inverse  problems  in  linear  isotropic  elasticity  are  considered.  The  problems 
of  interest  in  this  paper  are  Cauchy  problems  (see  [1])  that  cannot  be  solved 
directly;  therefore,  special  techniques  involving  transformations  are  required. 
The  regularization  method  employed  by  Gao  and  Mura  [2]  and  Yeih,  et.  al.  [3,  4] 
are  examined  in  order  to  improve  efficiency  and  computer  implementation.  At 
the  same  time,  an  alternative  method  based  on  the  classical  least-squares 
method  is  proposed. 

1.  Introduction 

Consider  the  following  problem.  A  thin  rectangular  plate  is  maintained  at  zero 
temperature  along  the  three  edges  a  =  0,  x  =  a,  and  v  =  0.  The  fourth  edge  y  =  is 


I 


subjected  to  a  temperature  distribution  f{x)  (see  Fig.  1).  We  would  like  to 
calculate  the  temperature  distribution  when  steady-state  is  reached. 

This  is  a  typical  Laplace  problem.  The  temperature  distribution  is  governed 
by  the  Laplace  equation, 


vv  = 


d^(p  d^(p 
dx^  dy^ 


=  0. 


(1) 


and  the  solution  is 


nnx  sinh(n;r>’/fl) 
a  sinh^rntbla)' 


(2) 


where 


2  .  .  .  nnx  . 

c,  =  - I  /(xjsin - dx, 

a  Jo  a 


(3) 


assuming  this  series  is  convergent. 

Now,  consider  the  following  problem.  Instead  of  prescribing  temperatures 
along  all  edges,  we  prescribe  a  temperature  f{x)  and  a  flux  along 

y  =  t>,  and  keep  the  remaining  edges  free  (see  [5])  as  shown  in  Fig.  2.  We  would 
like  to  calculate  the  temperature  distribution. 

Unlike  the  first  boundary-value  problem,  the  second  one  is  categorized  as  a 
Cauchy  (initial-value)  problem  (see  [1]),  and  it  is  a  type  of  inverse  problem.  The 
given  data  (boundary  conditions)  in  this  problem  are  called  the  Cauchy  data. 
Tikhonov,  one  of  the  first  pioneers  of  the  inverse  problem,  stated  tb?t  to  solve  an 
inverse  problem  means  to  discover  the  cause  of  a  known  result.  Therefore,  he 
continued  that  all  problems  of  the  interpretation  of  observed  data  are  actually 
inverse  [6]. 
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Although  his  description  seems  to  indicate  that  both  first  and  second 
problems  are  inverse;  however,  in  almost  all  cases,  we  distinguish  the  first  one 
from  the  second  one  and  call  the  first  one  a  forward  (regular)  problem  and  the 
second  one  an  inverse  problem.  The  inverse  problem  usually  involves  drawing  of 
a  conclusion  (cause)  from  incomplete  information.  Only  incomplete  information 
may  be  obtained  because  the  points  of  interest  are  visually  hidden,  or  the  sensors 
cannot  simply  reach  those  points. 

Two  of  the  first  classical  works  on  the  inverse  problem  may  be  attributed  to 
Legendre  [7]  and  Gauss  [8]  when  they  independently  calculated  the  trajectory  of 
Halley’s  comet  from  discrete  observed  data  by  means  of  the  least-squares 
method.  Today,  the  inverse  problem  is  receiving  more  attention,  particularly  in 
the  field  of  design  optimization  (e.g.,  Sobieczky  [9]  and  Staniz  [10]). 


2.  Inverse  Problems  in  Linear  Elasticity 

Let  us  see  how  the  inverse  problem  fits  in  linear  elasticity.  Historically,  inverse 
problems  in  elasticity  date  back  as  far  as  1907  to  Alamansi  [11].  A  more  recent 
work  was  done  by  Gao  and  Mura  [2]  and  Yeih,  et.  al.  [3,  4]. 

Consider  the  following  problem.  An  thin  circular  disc  of  a  radius  r  =  1.0m  and 
elasticity  constants  £  =  500MPa  and  v  =  0.3  is  subjected  to  a  radial  load 
/7  =  2.0MPa  as  shown  in  Fig.  3.  If  we  eliminate  rigid  body  motion,  the  x  and  y 
displacements  along  the  edge  are 


u,  = - cos(y), 

'250  '  ' 

w,  =  -^^sin(0). 
^  250  ^  ^ 


(4) 


3 


Now,  consider  the  following  problem.  The  same  disc  is  subjected  to  the  same 
radial  load  and  displacements  given  by  (4)  along  the  edge  ranging  0<6<n  as 
shown  in  Fig.  4.  Can  we  find  the  tractions  and  displacements  along  the 
remaining  edge  n<d<2n  assuming  that  the  disc  is  in  static  equilibrium? 
Obviously,  the  solutions  are 


p,  =2cos(e), 
=2sin{e). 


w,  = - cos(0), 

‘250  ^  ^ 

M,  = - sin(0). 

^250  ^  ^ 


(5) 


This  inverse  problem  is  simple  because  of  the  problem’s  simple  geometry  and 
loading  condition;  however,  how  can  we  treat  inverse  problems  that  have 
complicated  shapes  and  loading  conditions? 


3.  Numerical  Approach 

When  a  body  has  a  complicated  geometry  and  is  subjected  to  complex 
displacement  and  loading  conditions,  there  is  very  little  hope  to  solve  the 
problem  analytically.  In  general,  such  a  problem  is  solved  numerically  by  means 
of  the  finite  element  method  (FEM)  or  the  boundary  element  method  (BEM): 

jB^DBud^2  =  jNWa  FEM 

.  ^  "  (6) 

-Nu  + jTNui^r  =  JuNtrfr.  BEM 

^  r  r 
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normal  derivative,  Q  is  the  body,  and  F  is  its  boundary  which  is  assumed  to  be 
smooth  in  this  case  (the  presence  of 'A  coefficient).  Two  major  differences  between 
these  methods  is  that  FEM  solves  a  domain  integral  equation  whereas  BEM 
solves  a  boundary  integral  equation,  and  FEM  solves  for  displacements  only 
whereas  BEM  solves  for  both  displacements  and  tractions.  Both  methods  employ 
discretization  of  the  domain  of  integration;  consequently,  FEM  requires  to 
discretize  the  entire  body  whereas  BEM  requires  to  discretize  its  boundary  only 
when  body  forces  are  neglected. 

At  the  first  glance,  because  FEM  requires  prescribing  of  the  conditions  inside 
the  body,  this  method  is  not  suitable  for  our  purpose.  On  the  other  hand,  because 
BEM  requires  prescribing  of  the  conditions  along  the  boundary  only,  we  will 
employ  this  method  to  solve  Cauchy  problems. 

4.  The  General  Numerical  Procedure  to  Solve  Integral  Equations 

Before  proceeding  to  solve  Cauchy  problems,  we  first  look  into  the  procedure  to 

solve  an  integral  equation  numerically.  Consider  the  following  integral  equation, 
j^x'xu{x)dx-^u{x')  =  cos{nx'),  0<x'<l.  (7) 

The  exact  solution  is 

-  2  cos(  jrx\  (8) 

K 

and  we  would  like  to  approximate  this  solution  numerically. 

We  use  the  3-pt  Gauss-Legendre  quadrature  to  evaluate  the  integral  in  (7). 
The  formula  is 
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We  use  the  3-pt  Gauss-Legendre  quadrature  to  evaluate  the  integral  in  (7). 
The  formula  is 

fVWdx  =  y5^/(c  +  /^,  )w,,  (9) 

“  i=l 

where  J  =  {b-a)ll, 

c  =  (i>  +  a)/2, 

^  =  0,  V^},  (10) 

w  =  {5/9  8/9  5/9}. 

We  do  not  know  the  values  of  u  at  the  points  of  evaluation  (sample  points). 
Therefore,  we  have  three  unknowns  and  require  three  equations  to  solve  for 
them.  We  choose  the  same  points  for  x'  in  order  to  create  three  equations.  In 
essence,  we  have  transformed  an  integral  equation  into  a  system  of 
simultaneous  equations, 

‘  0.0278  -0.4965  0.0250][m,]  [0.9380 

0.1232  0.1565  -0.3889  J  0  >.  (11) 

-0.2813  1.1517  0.1972J[m3J  [-0.9380 

The  solution  is 

w  =  {-1.6030  1.2101,  4.0248}  (12) 

which  is  in  good  agreement  with  the  exact  solution, 

M  =  {-1.6019,  1.2159,  4.0336}.  (13) 

As  the  number  of  points  increases,  the  accuracy  of  solution  increases  as 
shown  in  Fig.  5.  The  rate  of  convergence  depends  on  the  numerical  integration 
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scheme.  Fig.  5  uses  the  trapezoidal  rule,  which  is  not  a  powerful  scheme,  in  order 
to  illustrate  this  point.  This  method  of  solving  integral  equations  is  called  the 
Nystrom  method  (see  [12]). 

5.  The  Compatibility  Condition  of  the  Cauchy  Data 

In  Section  4,  we  saw  how  an  ordinary  integral  equation  can  be  solved  by  the 
Nystrom  method.  Now  we  look  into  the  Cauchy  problem  and  pay  special 
attention  to  Laplace-type  problems  as  shown  in  Figs.  1  and  2.  But,  first,  we  list 
here  3  types  of  boundary- value  problem  (see  [1,  13]): 

i.  Dirichlet  Problem:  Only  (p  is  prescribed  along  the  boundary.  There  exists  a 
unique  solution. 

it  Neumann  Problem:  Only  the  normal  derivative  d(pldn  is  prescribed  along 
the  boundary.  The  solution  may  not  exist,  or  even  if  it  does,  it  is  not  unique  by  an 
arbitrary  constant. 

Hi.  Mixed  Problem:  Either  (p  or  d(p/dn  are  prescribed  along  the  boundary. 
There  exists  a  unique  solution. 

Bearing  in  mind  the  above  problem  types,  we  recall  Green’s  second  identity 
(see  [1]), 

(pix'}  =  j  ^  dr{x)  (14) 

where  G  is  the  Gieen’s  function  for  the  Laplace  equation.  This  identity  states 
that  when  both  cp  and  d(pldn  are  prescribed  along  the  boundary,  (p  inside  the 
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domain  can  be  uniquely  determined.  However,  this  statement  is  in  contradiction 
to  the  imiqueness  of  solution  of  the  Dirichlet  problem.  If  we  already  know  (p 
everywhere  along  the  boundary,  cannot  we  determine  (p  uniquely  without 
necessitating  the  knowledge  of  d(p/dn?  For  this  reason,  we  cannot  arbitrarily 
prescribe  (p  and  d<pldn  on  the  same  boundary  (Cauchy  data)  and  solve  the 

Cauchy  problem  directly.  This  difficulty  raises  the  compatibility  condition  which 
must  be  imposed  on  (p  and  d<pldn. 

Let  us  assume  that  the  boundary  F,  (p,  and  d(pldn  are  parametrically 
representable  by  the  variable  5,  so  that  we  can  write  r{s),(pis)  and  d(p{s)/dn.  In 
two  dimension,  the  normal  derivative  d<pis)l3n  is  written  as 


where 


„  -^y 

dx 

ds' 


(15) 


(16) 


Also,  we  can  write 


d(p  _d(p  dx  ^  d(p  dy 
ds  dx  ds  dy  ds 


(17) 


Therefore,  the  only  unknowns  are  d(p/dx  and  d(p/dy;  consequently  we  can  form  a 
system  of  simultaneous  equations. 


dx 

dy 

id<p] 

d<p 

ds 

dx 

ds 

dy 

dx 

d<p 

> = • 

dcp 

ds 

ds 

dn 

(18) 
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The  above  system  is  solvable  as  long  as  the  coefficient  matrix  is  nonsingular. 


We  then  proceed  to  compute  the  second  derivatives, 

d^(p  _  d  d(p 
dx^  dx  dx ' 

d^(p  d  d(p 
dy^  dy  dy 

We  take  advantage  of  the  chain  rule  again  to  form, 


^  ' 
ds  dx 

]A 

ds  [(9y 


whence 


\  ^  ] 

dx 

d 

j 

dy  d 
ds  ds 
dx 

ds  dn 


where 


[ds  [ds 


The  computed  second  derivatives  must  satisfy  the  Laplace  equation. 


I  -n 

dx^  dy^  ■ 


The  Cauchy  data  are  then  compatible. 
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6.  Numerical  Difficulties  in  Solving  Cauchy  Problems 

Consider  a  domain  Q  which  is  bounded  by  a  boundary  F.  The  Laplace  equation 
holds  inside  this  domain.  Denote  a  part  of  the  boundary  as  7^  and  the  remaining 
part  as  Fj,  so  that  F=rj  +  7^.  Suppose  the  Cauchy  data  (potential  and  flux)  are 
prescribed  along  We  would  like  to  compute  the  potential  and  flux  along 

For  the  simplicity  sake,  we  assume  that  the  Cauchy  data  are  compatible,  and 
with  the  Cauchy  data,  we  can  carry  out  the  integrations  along  F^  in  (14)  and 
denote  the  sum  as  b{x'), 

b{x')  =  f  -  .^3^  <p{x)\dF.  (24) 

•  on  on 

1,1-  -I 

Consequently,  Green’s  second  identity  is  reduced  to 

vW)  +  <pix)  -  a(x.x')^^^dr = b(x’),  x’^n.  (26) 

r  dn  on 

/jU  j 

We  do  not  wish  to  keep  x'  inside  Q  since  we  do  not  have  any  information  as  to 
the  value  of  (p.  When  x'  is  placed  on  the  boundary  F  and  if  the  boundary  is 
smooth,  (25)  can  be  rewritten  as 

1  +  f _  (;(;c,x')^^]dF  =  b{x\  x'  e  F.  (26) 


We  apply  the  Nystrom  method  here.  Let  us  use  the  Nth  order  numerical 
integration  scheme  to  rewrite  (26)  as 


e  F,  i  =  1,. 
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From  (27),  we  see  that  we  can  construct  one  equation  at  each  sample  point. 
However,  there  are  2  unknown  quantities,  (p  and  d(p/dn  at  each  sample  point,  so 

that  there  are  more  unknown  than  the  number  of  equations.  In  general,  we 
cannot  solve  this  system  of  equation  directly.  The  dimensions  of  the  coefficient 
matrix  are  N  by  2N,  and  even  if  all  the  equations  are  independent,  there  is  an 
2N  -N  =  N  different  dimensional  family  of  solutions.  There  is  no  known  method 
to  determine  the  correct  solution  directly. 

There  is  a  way  to  reduce  the  order  of  dimensional  family  of  solutions,  however. 
Since  there  is  nothing  to  prevent  us  from  placing  x'  in  (27)  on  fj,  we  can  increase 

the  number  of  equations  without  increasing  the  number  of  unknowns.  In  the 
present  case,  if  we  choose  N  points  on  7^,  we  can  have  an  2N  x2N  system  of 

equations  which  we  may  be  able  to  solve  directly.  It  should  be  pointed  out  here 
that  although  we  seem  to  construct  a  well-balanced  system  of  equations  by  this 
scheme,  it  is  not  always  the  case  as  we  will  see  later. 

7.  The  Ref^larization  Method 

When  a  system  of  simultaneous  equations  is  ill-conditioned,  that  is,  the 
coefficient  matrix  is  nonsquare  or  singular,  we  cannot  solve  it  by  means  of 
Gaussian  elimination  for  example.  In  general,  we  cannot  expect  an  extremely 
good  solution  from  an  ill-conditioned  system;  however,  we  can  obtain  a  probable 
solution  by  imposing  extra  conditions  as  to  the  nature  of  the  solution.  Gao  and 
Mura  [2]  introduced  the  regularization  method  by  applying  an  integral 
transformation  to  the  original  integral  equation  (e.g.,  (26)).  Although  their 
transformation  is  theoretically  valid,  computer  implementation  of  the  resulting 
multiple  integral  equation  is  difficult  and  not  numerically  efficient.  It  is  more 
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convenient  if  we  apply  a  linear  algebraic  transformation  to  the  numerical 
equation  (e.g.,  (27)).  Let  us  write  (27)  in  the  linear  algebraic  form, 

Ax  =  b  (28) 

where  the  vector  x  contains  the  unknown  (p  and  d<p/dn.  Because  A  is  either 

nonsquare  or  singular,  its  inverse  does  not  exist.  The  approach  Gao  and  Mura  [2] 
took  was  to  minimize  the  norm  of  x  subject  to  (28)  or 

minimize  Ixjlj  (29) 

subject  to  ||Ax  -  bjjj  <  e,  (30) 

where  11*112  =  Vx^x.  (31) 

The  constant  e  is  user-defined.  The  Lagrangian  function  for  this  problem  is 

L(x;a)  =  ax^x-(Ax-b)’^(Ax-b).  (32) 

The  corresponding  Euler  equation  (see  [14])  is 

(A^A  al)x  =  A^'b,  (33) 

so  that  the  solution  is 

X  =  (A^A  +  al)  '  A^b.  (34) 

The  value  of  a  has  to  be  determined  recursively  until  (30)  is  fulfilled. 

The  computer  implementation  of  (34)  is  far  easier  than  that  of  Gao  and 
Mura’s  multiple  integral  equation.  However,  computing  A^A  is  susceptible  to 
numerical  roundoff  errors  and  not  generally  recommended.  If  A  is  a  full  rank 


12 


matrix,  the  QR  decomposition  (see  [15])  maj'  be  helpful,  but  we  would  like  to 
consider  a  more  general  case.  We  introduce  here  the  singular-value 
decomposition  (SVD)  technique  (see  [15]).  Basically,  the  SVD  decomposes  A  into 

A  =  Ul:v^  (35) 


where  the  elements  of  the  diagonal  matrix  Z  are  called  the  singular  values  of  A. 
The  matrices  U  and  V  are  unitary  matrices,  that  is. 


U^'U  =  I, 
\^\  =  l 


(36) 


After  substituting  (35)  into  (33),  we  obtain 

(Z+aZ-’')Vx  =  U^b, 


(37) 


and  the  solution  (34)  becomes 

x  =  V(Z+aZ'^)'‘u’'b.  (38) 

The  solution  by  (38)  is  easier  to  control  because  we  are  tampering  with  the 
singular  values  only  which  directly  affect  the  condition  number  of  A.  Also, 
although  two  matrix  inversion  are  involved  in  (38),  because  the  matrix  Z  is  a 
diagonal  matrix  so  that  it  is  a  matter  of  inverting  the  diagonal  elements.  The 
regularization  method  coupled  with  the  SVD  is  recommended  to  deal  with  all 
kinds  of  Cauchy  problem. 
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8.  The  Least-Squares  Method 

The  least-squares  method  has  been  used  for  many  inverse  and  optimization 
problems.  Legendre  [7]  and  Gauss  [8]  used  this  method  to  solve  one  of  the  first 
inverse  problems.  We  employ  this  method  here  coupled  with  the  SVD  introduced 
in  the  previous  section. 

Given  an  equation  Ax  =  b,  the  least-squares  method  is  derived  as  follows.  We 
would  like  to  minimize  (extremize) 

llAx-bllj  (39) 

with  respect  to  x.  Thi.s  gives  rise  to 

— (x^A^Ax -2x^A^b-i-b^b)  =  0,  (40) 

or 

A’‘Ax  =  A^b.  (41) 

This  equation  is  known  as  the  normal  equation.  Commonly,  (41)  is  solved  by 
means  of  the  QR  decomposition  (see  [15])  because  of  the  suceptability  of  the  A^A 
computation  to  roundoff  error.  The  decomposed  normal  equation  becomes 

Rx  =  Q^b  (42) 

The  above  equation  is  merely  a  backsubstitution  operation. 

The  use  of  the  QR  decomposition  is  valid  as  long  as  an  m  x  n  matrix  A  has  a 
rank  n,  so  that  A^A  is  a  full  rank  matrix.  However,  many  inverse  problems  have 
a  less  than  rank  n  matrix.  Consequently,  the  QR  decomposition  cannot  be  used. 
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In  these  cases,  we  employ  the  singular  value  decomposition  technique.  The  SVD 
decomposes  A  into 


A  =  VEL!'. 

where  V^V  =  I, 

U^U  =  I. 


Then  (41)  is  reduced  to 


ZV^x  =  U^'b. 


Consequently, 


X  =  VZ  U’^b. 


(43) 

(44) 


(45) 


(46) 


Notice  the  resemblance  between  (38)  and  (46).  If  we  let  a  in  (38)  go  to  zero,  the 
two  equations  are  the  same.  Therefore,  the  regularization  method  and  the 
least-squares  method  are  very  similar;  however,  they  are  not  the  same.  While  the 
regularization  method  attempts  to  balance  the  order  of  singular  values,  the 
least-squares  method  cuts  off  numerically  negligible  singular  values.  Therefore, 
the  regularization  method  is  a  full  expansion  whereas  the  least-squares  method 
is  a  tuncated  expansion. 


9.  Linear  Elasticity  Problems 

We  now  apply  the  least-squares  method  to  some  two  dimensional  isotropic  linear 
elasticity  problems.  The  basic  equation  we  use  here  is  known  as  Somigliana’s 
boundary  integral  equation  (see  [16]), 


(47) 


^  (x')  +  /  7;  {x,x')uj  (x)  ds{x)  -fu,j  {x,x')tj  ix)ds{x). 


where  Uij{x,x')  = 


-1 


8;r(l-v)G 


{(3-4v)ln(r)5^-r.rJ, 


"  4;Rr-v)r{t**  -  )}■ 


n  =  Xi-x„ 

__  Br  _r- 
dx^  r ' 

n,  =  Normal  vector, 

=  Kronecker  delta, 
G  -  Shear  modxilus, 

V  =  Poisson’s  ratio. 


(48) 


The  boundary  F  is  assumed  to  be  smooth  as  shown  in  Fig.  6  because  the  constant 
'A  depends  on  the  smoothness  of  the  boundary. 

Let  us  consider  a  mixed  boundary-value  problem  first.  For  the  sake  of 
simplicity,  let  us  suppose  that  displacements  are  prescribed  along  a  part  of  F 
denoted  by  and  tractions  are  prescribed  along  the  remaining  boundary 
denoted  by  Fj,  so  that  F=  +  f^.  After  we  move  the  unknown  displacement  and 
traction  terms  on  the  left  side  of  (47)  and  the  known  displacement  and  traction 
terms  on  the  right,  we  obtain 


]-u,{x')-^  \T^,{x,x%{x)ds{x)-  \u,^{x,x')t,{x)ds{x) 

^  r,  r, 

=  -jF^(x,x>^(x)d:s(x)+  ^U^^{x,x')t^{x)ds{x),  (49) 

r,  r. 
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where  the  position  of  the  first  term  on  the  left-hand  side  is  not  necessarily  on 
this  side.  If  x'  is  on  7^,  this  term  moves  over  to  the  right-hand  side. 

In  general,  these  integrations  are  carried  out  numerically,  and  the  one  of  the 
most  common  techniques  associated  with  integration  is  to  discretize  (subdivide) 
the  domain  of  integration  into  smaller  subdomains  or  elements.  Then  (49)  is 
rewritten  as 


1  **  t 


p(") 
*2 


Asl  j^m) 


M+N 


n=M+l  rl*5 

*  ■i 


ms]  p(«> 


in  which  7^  and  are  divided  into  M  and  N  elements,  respectively.  We  apply  a 
Qth  order  numerical  integration  technique  to  (50),  but  instead  of  giving  the 
values  of  the  displacements  and  the  tractions  at  each  sample  point  exactly,  we 
approximate  these  values  by  based  on  P  values  in  each  element,  that  is,  we 
appriximate  the  displacements  and  tractions  in  each  element  by 


m=l 


(51) 


where  the  function  5  is  called  the  shape  function.  After  we  substitute  (51)  into 
(50),  we  carry  out  the  integration  and  obtain 
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(52) 


n=M4*l  k-\  m=l 


In  the  case  of  the  constant  approximation  (P  =  l),  there  are  2{M  +  N)  unknowns 
(displacements  and  tractions)  as  shown  in  Fig.  7,  and  we  require  the  same 
number  of  equations  in  order  to  solve  for  them.  Because  at  each  x'  we  can 
construct  two  equations,  we  choose  x'  at  the  same  points  as  x  (see  Fig.  7). 
Moreover,  because  the  displacements  and  tractions  are  given  on  the  right-hand 
side  of  (52),  we  can  simplify  the  right-hand  side  terms.  If  we  write  (52)  in  matrix 
form  it  becomes 


(53) 


or  even  more  simply, 


where 


Ax  =  b, 

A=[T  -U], 
x  =  {u 


(54) 

(55) 


If  the  coefficient  matrix  A  is  well-conditioned,  x  can  be  solved  by  a  simple 
technique  like  the  Gauss  elimination.  This  is  the  starndard  boundary  element 
method  solution. 


10.  Cauchy  Problems  in  Linear  Elasticity 

We  now  examine  some  Cauchy  problems  in  isotropic  linear  elasticity.  Let  us  first 
consider  a  general  problem.  We  can  modify  the  problem  in  the  previous  section 
and  state  that  both  displacements  and  tractions  are  prescribed  along  and  no 
values  are  prescribed  along  7^  (free  boimdary).  Then  we  have  to  modify  (49)  as 

^uXx')-^  jT,jix,x')uj(x)ds{x)-  ju,jix,x')tj{x)dsix) 

=  -  J  T^{x,x')ujix)ds{x)  +  j  U,j{x,x')tjix)ds{x).  (56) 

r,  r, 

When  (56)  is  discretized  as  (49)  was  discretized  into  (52)  then  into  (54),  it  again 
becomes  the  same  form  as  (54).  However,  under  the  constant  assumption  the 
total  nximber  of  unknowns  in  this  case  is  4N  (2  displacements  and  2  tractions  at 
each  point),  and  the  total  nixmber  of  equations  is  still  2{M  +  N).  Therefore,  the 

coefficient  matrix  A  is  nonsquare,  and  we  cannot  employ  a  simple  technique  to 
solve  the  simultaneous  equations. 

We  might  be  tempted  to  choose  M  =  N  so  that  the  coefficient  matrix  to  be 
square,  but  this  scheme  usually  does  not  work  for  the  following  reason.  Suppose 
fj  covers  only  1%  of  the  entire  boundary  and  the  free  boundary  covers  the  rest. 
In  this  case,  even  if  we  construct  a  square  coefficient  matrix  by  taking  M  =  N,  it 
is  still  highly  singular  unless  the  elastic  domain  and  the  loading  have  a  high 
degree  of  symmetry.  In  general,  however,  this  is  not  the  case. 

To  overcome  the  difficulty  of  constructing  a  nonsingular  square  matrix  we 
require  special  techniques  such  as  the  regularization  method  and  the 
least-squares  method.  Once  the  integral  equation  (56)  has  been  converted  into 
the  form  (54),  the  procedure  is  identical  as  explained  in  Sections  7  and  8. 
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Sample  Problem  1 

We  employ  the  least-squares  method  discussed  in  Section  8  to  solve  the  following 
problem.  The  upper  half  of  a  circular  disk  of  radius  Im  is  imder  radial  pressure 
of  2MPa  as  shown  in  Fig.  4.  The  center  of  the  disk  is  fixed  and  its  Young’s 
modulus  and  Poisson’s  ratio  are  500MPa  and  0.3,  respectively.  The 
displacements  along  the  upper  half  boundary  are  given  as 


“x  =^cos(0)m, 

Uy  =  (0  <  0  <  7t). 


(58) 


From  these  information,  we  would  like  to  compute  the  displacements  and 
tractions  along  the  boundary  of  the  lower  half  of  the  disk. 

The  exact  solutions  are 


“x  =^cos(0)m, 

=3^sin(0)m,  (59) 

=2  cos  (5)  MPa, 

ty  =  2sin(ff)  MPa,  (^<0<  2n). 

We  use  40  equal  length  constant  elements  for  the  prescribed  and  free 
boundaries.  The  prescribed  boundary  is  along  (0  <  ^  <  ;r)  and  the  free  boundary 
r ,  is  along  (n<d<  2n).  After  we  substitute  them  into  (56)  and  convert  into  a  form 
of  (54),  and  apply  the  least-squares  method,  we  obtain  the  results  as  shown  in 
Figs.  8-11. 

We  notice  that  even  though  the  solutions  seem  to  agree  with  the  exact 
solutions,  there  are  still  significant  discrepancies  between  them.  We  investigate 


20 


cause  of  this  discrepancies  a  little  further.  Suppose  we  solve  a  mixed 
boundary-value  problem  of  the  same  elastic  body  with  the  following  boundary 
conditions: 

M,  =  0  at  0  =  njl  and  3«/2, 

M,  =  0  at  0  =  0  and  n, 

(60) 

r,  =2  cos  (6)  MPa, 

ty  =  2 sin  (6)  MPa,  (0  <  6  <  2;r). 

We  use  the  same  discretization  as  before.  Then  we  obtain  the  results  as  shown  in 
Fig.  12.  It  seems  to  indicate  that  for  the  current  discretization  the  displacements 
in  (58)  are  not  compatible  with  the  given  tractions.  It  seems  to  make  more  sense 
to  prescribe  the  displacements  obtained  from  this  mixed  boundary-value  problem 
in  place  of  (58).  The  results  of  solving  this  Cauchy  problem  are  shown  in  Figs. 
13-16.  The  restilts  are  in  better  agreement. 

Sample  Problem  2 

A  side  (x  =  1/2)  of  a  thin  square  plate  of  side  Im  is  pulled  by  a  uniform  traction  of 
10,000Pa  as  shown  in  Fig.  17.  The  displacements  along  this  side  is  measured 
Yoimg’s  modulus  and  Poisson’s  ratio  are  30MPa  and  0.25,  respectively.  We  would 
like  to  compute  the  displacements  and  tractions  along  the  remaining  sides. 

In  this  problem,  the  measured  displacements  are  the  displacement  solutions 
for  the  mixed  boundary-value  problem  in  which  the  boundary  conditions  are 
given  by 
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(61) 


w,  =  0  at  jr  =  0, 

=  0  everjrwhere 

t,  =100Paat(;c  =  ±V2,-l/2<y<l/2), 
r,  =  0  at  (-1/2  <  x  <  1/2,  y  =  ±  1/2), 
r,  =0at(;c  =  ±V2,-l/2<y<l/2). 

Unlike  the  circxilar  disk  in  the  previous  problem,  this  plate  has  sharp  comers 
(geometric  singularities).  We  would  like  to  see  how  the  least  squares  method 
works  with  this  type  of  irregularity  (discontinuity).  Figs.  18-21  show  the  results, 
and  we  see  that  the  results  are  not  as  good  as  ones  for  the  previous  problem.  The 
farther  away  is  the  point  of  interest,  the  worse  become  the  solutions. 

Sample  Problem  3 

A  circular  elastic  disc  embedded  in  another  thin  elastic  body  and  loaded  as 
shown  in  Fig.  22.  The  boiindary  of  the  outer  elastic  body  is  composed  of  straight 
lines  denoted  by  ®,  (D,  ®,  and  ®.  Suppose  the  Cauchy  data  (displacements  and 

tractions)  are  prescribed  along  ®  and  (D.  We  would  like  to  compute  the  remaining 
boundary  values.  The  Cauchy  data  are  the  displacements  and  tractions 
computed  by  solving  the  corresponding  mixed  boundary-value  problem.  A  total  of 
144  constant  boundary  elements  are  used.  The  resulting  deformation  is  shown  in 
Fig.  23.  The  Cauchy  solution  agree  with  the  mixed  boimdary- value  solution. 

Conclusions 

Cauchy  problems  in  linear  isotropic  elasticity  are  solvable  although  not  directly 
by  means  of  the  regularization  method  as  well  as  the  least-squares  method 
coupled  with  the  singular  value  decomposition.  Unlike  the  regular 
boundary-value  problem,  the  Cauchy  problem  is  very  sensitive  to  compatibility 


among  the  Cauchy  data.  Feeding  the  exact  solution  does  not  necessarily 
guarantee  a  good  solution.  The  accuracy  depends  heavily  on  how  the  body  in 
equation  is  discretized.  The  Nystrom  type  solution  simplifies  implementation  on 
computer  and  improves  performance. 
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Figure  Captions 

Figure  1.  A  regular  boundary- value  problem. 

Figure  2.  A  Cauchy  problem. 

Figure  3.  A  regular  elasticity  boundary-value  problem. 

Figure  4.  A  elasticity  Cauchy  problem. 

Figure  5.  The  Nystrom  solution  comparison. 

Figure  6.  The  parameters  in  Somigliana’s  boundary  integral  equation. 
Figure  7.  A  discretized  boundary. 

Figure  8.  The  Cauchy  solution  for  the  x  displacement. 

Figure  9.  The  Cauchy  solution  for  the  y  displacement. 

Figure  10.  The  Cauchy  solution  for  the  x  traction. 

Figure  11.  The  Cauchy  solution  for  they  traction. 

Figure  12.  The  inherent  errors  in  discretization. 

Figure  13.  A  more  compatible  Cauchy  solution  for  the  x  displacement. 
Figure  14.  A  more  compatible  Cauch;  solution  for  the  y  displacement. 
Figure  15.  A  more  compatible  Cauchy  solution  for  the  x  traction. 

Figure  16.  A  more  compatible  Cauchy  solution  for  the  x  traction. 

Figure  17.  A  thin  square  plate  under  load. 

Figure  18.  The  Cauchy  solution  for  the  x  displacement. 

Figure  19.  The  Cauchy  solution  for  the  y  displacement. 

Figure  20.  The  Cauchy  solution  for  the  x  traction. 

Figure  21.  The  Cauchy  solution  for  they  traction. 

Figure  22.  A  more  complicated  problem. 

Figure  23.  The  Cauchy  solution  for  deformation. 
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Load  Diffusion  and  Absorption  Problems  from  a  Finite  Fiber  to  Elastic  Infinite 
Matrix  by  the  Equivalent  Inclusion  Method 

Ven-Gen  Lee  and  Toshio  Mura 
Theoretical  and  Applied  Mechanics 
Northwestern  University 
Evanston,  III  60208 

Abstract 

The  load  transfer  behavior  of  a  finite  fiber  perfectly  bonded  to  an  infinite  matrix  of 
distinct  elastic  properties  is  investigated  in  this  paper.  The  fiber  is  subjected  to  the 
uniform  distributed  loading  applied  at  infinity  or  on  one  cross  section  of  the  fiber.  The 
stress  disturbance  due  to  the  existing  fiber  is  simulated  by  the  equivalent  inclusion 
method,  which  formulates  the  inhomogeneity  problem  to  a  system  of  integral  equations. 
By  dividing  the  fiber  into  finite  numbers  of  ring  elements  with  uniform  distributed 
eigenstrains,  the  integral  equations  can  be  further  reduced  to  a  system  of  algebraic 
equations  with  coefficients  expressed  in  terms  of  the  integrals  of  Lipschitz-Hankel  type. 
Numerical  results  are  presented  for  resultant  axial  force  for  various  fiber  length  and 
material  properties.  The  limiting  cases  of  the  infinite  and  semi-infinite  fibers  are  also 
compared  with  the  exact  and  approximate  solutions. 
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